Numerical mode terrain processing method based on various fusion data

Abstract

The invention discloses a numerical mode terrain processing method based on various fusion data, and relates to the technical field of geographic space industrial data processing, and the method comprises the steps: calculating a gradient, a curvature and a landform boundary line based on the preprocessed multi-source terrain grid data, forming a structure guide field, and carrying out the construction of the structure guide field; establishing a confidence coefficient model to generate a confidence coefficient field based on the sensor error, the spatial density and the landform consistency score, and generating a unified landform field by adopting a structure-preserving fusion algorithm; according to the data difference between the unified terrain field and the preprocessed multi-source terrain grid data and in combination with the geometric continuity of the terrain, identifying a local offset region with abnormal gradient, gradient abrupt change and boundary mismatch, and calculating an offset error field by adopting local registration and gradient consistency constraint; and error compensation and local reconstruction are carried out on the unified terrain field to obtain a corrected terrain field. The numerical mode initial field quality and the stability and usability of the simulation result are improved.

Description

Technical Field

[0001] This invention relates to the field of geospatial industrial data processing technology, and in particular to a numerical model terrain processing method based on multiple fused data. Background Technology

[0002] Topographic data, as a crucial input for the initial field and underlying surface parameterization of numerical models, is widely used in scenarios such as numerical weather prediction, climate simulation, hydrological runoff and flood projection, pollution diffusion assessment, and wind energy resource calculation. With the increasing demand for high-resolution models, refined urban simulations, and environmental protection in industrial parks, topographic data acquisition has gradually expanded from a single digital elevation model (DEM) to a multi-source parallel system, including various data formats such as DEMs, laser point clouds, remote sensing inversion topography, and surveyed vector topography. For engineering applications, the relevant processes typically organize topographic data flow using industrial data processing methods. This involves registering and managing metadata for multi-source data, unifying coordinate and vertical reference standards, performing noise suppression, anomaly detection, and missing data processing, followed by rasterization and grid alignment. Structural features such as slope, curvature, and boundary response are then extracted to support subsequent fusion, correction, and grid adaptation, thereby forming a topographic field that meets the input format and operational constraints of numerical models.

[0003] Existing technologies still have shortcomings in multi-source fusion and model grid adaptation, mainly in two aspects: First, multi-source data differ significantly in resolution, coverage density, error structure, and boundary representation. Although quality information can be introduced for weight allocation and structural constraints can be applied during fusion processing, problems such as weakened structural features, boundary drift, local morphological distortion, and loss of geometric continuity can still easily occur in complex terrain boundaries and local detail areas. Second, when resampling and smoothing different numerical model grid structures such as regular grids, hexagonal grids, and nested grids, it is difficult to maintain a stable balance between gradient control in steep regions and global smoothing. This can easily lead to unstable initial field gradients and induce integral noise, affecting the stability and usability of simulation results. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a numerical model terrain processing method based on multiple fused data to solve the problems of difficulty in balancing structural consistency and geometric continuity in the multi-source fusion stage, and difficulty in maintaining stable gradient control and global smooth balance in the model grid adaptation stage.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution: This invention provides a numerical model terrain processing method based on multi-source fusion data. The method includes: acquiring multi-source terrain data, preprocessing it, and mapping it to a unified reference grid to generate preprocessed multi-source terrain grid data; calculating slope, curvature, and landform boundary lines based on the preprocessed multi-source terrain grid data to form a structure-guided field; establishing a confidence model based on sensor error, spatial density, and landform consistency score to generate a confidence field; and using a structure-preserving fusion algorithm to generate a unified terrain field; identifying local offset regions with slope anomalies, gradient abrupt changes, and boundary mismatches based on the data differences between the unified terrain field and the preprocessed multi-source terrain grid data, combined with terrain geometric continuity; calculating the offset error field using local registration and gradient consistency constraints; and performing error compensation and local reconstruction on the unified terrain field to obtain a corrected terrain field; performing model grid adaptation and terrain resampling based on structure-guided field constraints on the corrected terrain field according to the numerical model grid structure to generate a grid terrain field; and obtaining a smoothed terrain field by applying a slope constraint strategy to steep areas within the grid terrain field and then performing overall smoothing processing.

[0007] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the multi-source terrain data includes digital elevation model data, laser point cloud data, remote sensing inversion terrain data, and surveying vector terrain data.

[0008] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for generating preprocessed multi-source terrain grid data are as follows: According to the data registration records, the multi-source terrain data is subjected to unified coordinate benchmark transformation and elevation benchmark conversion, and noise identification and filtering are performed to generate preprocessed multi-source terrain data; The preprocessed multi-source terrain data is mapped to a unified reference grid to generate source elevation grids and overlay marker gratings. Terrain gap filling and outlier detection are performed on the source elevation grid to generate preprocessed multi-source terrain grid data.

[0009] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for forming the structure-guided field are as follows: Extract source elevation grids, overlay marker gratings, missing measurement masks, and filling masks from the preprocessed multi-source terrain grid data, and establish an effective cell set; Within the effective cell set, slope and curvature are calculated to generate slope and curvature rasters. Boundary response rasters are constructed on a unified reference grid to extract the terrain boundary line mask. The slope grid, curvature grid, and terrain boundary mask are combined into a structure guidance field using a unified reference grid index.

[0010] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for establishing a confidence model based on sensor error, spatial density, and terrain consistency score are as follows: An error baseline quantity is constructed on a unified reference grid and converted into an error confidence level through a monotonically decreasing mapping, thereby generating an error confidence level layer. Density baseline quantities are calculated on a unified reference grid and converted into density confidence values ​​through a monotonically increasing mapping, generating a density confidence layer. Based on the structure-guided field, the morphological difference quantity is calculated on a unified reference grid, and the boundary coincidence mark is calculated on the grid cell covered by the landform boundary line mask. The morphological difference quantity and the boundary coincidence mark are combined into a consistency benchmark quantity and the consistency confidence is obtained through difference reduction mapping, and a consistency confidence layer is generated.

[0011] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for generating a unified terrain field are as follows: The error confidence layer, density confidence layer, and consistency confidence layer are weighted and summed according to the combination coefficients to generate a comprehensive confidence score, which is then normalized to obtain the confidence field. The source elevation grid is weighted and synthesized based on the confidence field to obtain an initial unified topographic field. A structure-preserving fusion algorithm is then used to perform fusion iterations to generate a unified topographic field.

[0012] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for identifying local offset regions with slope anomalies, gradient abrupt changes, and boundary mismatches are as follows: Within the grid cells marked as valid by the overlay marker raster, calculate the difference raster between the source elevation grid and the uniform topographic field, and aggregate them into a comprehensive difference raster; The baseline slope raster and baseline gradient raster are calculated based on the unified topographic field, and combined with the comprehensive difference raster, slope anomaly indicators, gradient abrupt change indicators and boundary mismatch indicators are identified to generate a set of offset risk layers. Thresholding and connected component aggregation are performed on the set of offset risk layers to generate a local offset region mask.

[0013] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for obtaining the corrected terrain field by performing error compensation and local reconstruction on the unified terrain field are as follows: Local registration is performed within the area covered by the local offset region mask to obtain a set of local displacement vectors, and the set of local displacement vectors is interpolated to generate the initial displacement field; The displacement field is updated by applying gradient consistency constraints to generate the offset error field. Error compensation is performed on the uniform topographic field based on the migration error field, and local reconstruction is performed in the local migration area to generate the corrected topographic field.

[0014] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for generating the grid terrain field are as follows: Collect numerical model grid configuration parameters, establish numerical model grid structure, and generate a set of model grid indexes; Based on the corrected terrain field, terrain resampling is performed on the pattern grid index set to obtain the grid terrain field, and the structure guiding field is used as a constraint to obtain the grid terrain field that fits the pattern grid.

[0015] As a preferred embodiment of the numerical model terrain processing method based on multiple fused data described in this invention, the specific steps for obtaining the overall smoothed terrain field are as follows: Calculate the slope grid based on the grid terrain field, generate a steep region mask, and perform connected component aggregation to obtain a set of steep regions; Apply a slope constraint strategy within a set of steep regions to generate a slope-constrained grid terrain field; Based on the slope-constrained grid topographic field, the global grid is smoothly updated with the goal of stabilizing the initial gradient of the model, resulting in an overall smoothed topographic field.

[0016] The beneficial effects of this invention are as follows: by forming a set of numerical model terrain processing workflows for geospatial industrial data processing, multi-source terrain data can be fused and corrected under the premise of structural consistency, continuity and stability, and model grid adaptation can be completed. This makes the output terrain field more consistent in terms of landform boundaries, morphological details and gradient stability, and has a lower risk of local offset and mismatch, thereby improving the quality of the initial field of the numerical model and the stability and usability of the simulation results. Attached Figure Description

[0017] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 This is a flowchart of a numerical model terrain processing method based on multiple fused data.

[0019] Figure 2 The confidence level binning calibration curve is shown.

[0020] Figure 3 This is a curve showing the monitoring data during iterative convergence.

[0021] Figure 4 This is a graph showing the cumulative distribution of absolute error across the entire domain. Detailed Implementation

[0022] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0023] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0024] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0025] Reference Figures 1-4 This is one embodiment of the present invention, which provides a numerical model terrain processing method based on multiple fused data, including the following steps: S1: Collect multi-source terrain data, preprocess it, and then map it to a unified reference grid to generate preprocessed multi-source terrain grid data; S1.1: Collect multi-source terrain data and record the spatial reference system, vertical datum, resolution, acquisition time and coverage, and generate data registration records; Furthermore, digital elevation model data, laser point cloud data, remote sensing inversion terrain data, and surveying vector terrain data are collected, and a data source identifier is written for each data source; data header information and accompanying metadata are read, and the coordinate system name, projection type, and projection parameters are written into the spatial reference system field; the elevation datum type, datum surface parameters, elevation unit, and positive elevation direction are written into the vertical datum field; raster cell size, point density, average point spacing, effective cell ratio, contour interval, and scale level are written into the resolution field; and the collection start and end time and version time are written into the acquisition time field; the data boundary point sequence is extracted and converted to the coordinate system indicated by the spatial reference system field; the minimum bounding rectangle algorithm is used to calculate the four boundaries and write them into the coverage field; the target plane coordinate system, target vertical datum, and unified reference grid configuration in the numerical model grid configuration parameters are read and written into the same data registration record; checks are performed on missing fields, field format, and unit consistency, and the check results are written into the verification mark.

[0026] S1.2: Perform unified coordinate benchmark transformation and elevation benchmark conversion on the multi-source terrain data according to the data registration records, and perform noise identification and filtering to generate preprocessed multi-source terrain data; Furthermore, the spatial reference system field and the target plane coordinate system are read according to the data registration records. Projection forward and inverse calculations are used to convert the plane coordinates of each data source to the target plane coordinate system, and the coordinate axis units are written back according to the target plane coordinate system. The vertical datum field and the target vertical datum are read according to the data registration records. The elevation of each data source is converted to the target vertical datum using the datum surface difference table interpolation, and the elevation units and positive direction are written back according to the target vertical datum. Noise markers are written on the converted data and filtered out. For digital elevation model data, neighborhood median filtering is used to superimpose slope consistency check markers for spikes and burrs, and the markers are replaced with neighborhood medians. For laser point cloud data, stepwise morphological filtering is used to extract ground points and outliers are removed according to statistical outlier criteria. For remote sensing inversion topographic data, quality marker masks are used to remove low-confidence pixels. For surveying vector topographic data, topological consistency checks are used to remove duplicate and fault elements.

[0027] S1.3: Map the preprocessed multi-source terrain data to a unified reference grid to generate source elevation grids and overlay marker gratings; Furthermore, the unified reference grid configuration in the data registration record is read to establish a spatial range, grid step size, and grid index system. The preprocessed digital elevation model data is written into the grid cells using bilinear resampling. The laser point cloud ground points are aggregated by grid cells and written with median elevation. The effective pixels of the remote sensing inverted terrain are aggregated by grid cells and written with median elevation. The surveyed vector terrain features are rasterized and written into the grid cell elevation values. When there are multiple surveyed vector terrain features in the same grid cell, the writing source is selected according to the priority of scale level and the order of acquisition time, and a source mark is written. For each grid cell, a coverage mark raster and a noise mark raster are written simultaneously. A coverage mark value of 1 indicates that the cell has been written with elevation values, and a coverage mark value of 0 indicates that the cell has not been written with elevation values. A noise mark value of 0 indicates that the elevation value written to the cell comes from the data retained after filtering in S1.2, and a noise mark value of 1 indicates that the data corresponding to the cell has been hit by noise discrimination and has been removed or replaced. When the coverage mark value is 0, the noise mark is written as an invalid mark and is not counted in subsequent statistics.

[0028] S1.4: Perform terrain gap filling and outlier detection on the source elevation grid to generate preprocessed multi-source terrain grid data; Furthermore, a missing measurement mask is written based on the coverage marker raster, with a missing measurement mask value of 1 corresponding to a coverage marker value of 0. Outlier markers are written within grid cells where the coverage marker value is 1 using the Hamper filter rule, with the Hamper filter constructing outlier discrimination based on the neighborhood median and median absolute deviation. For grid cells where the outlier marker value is 1, inverse distance weighted interpolation is used to write back the elevation value, and the coverage marker for that grid cell is kept at 1. For grid cells where the missing measurement mask value is 1, inverse distance weighted interpolation is used to write back the filled elevation and write it into the filled mask, while simultaneously writing back the coverage marker for that grid cell to 1. The source elevation grid, coverage marker raster, missing measurement mask, filled mask, and outlier markers are then aggregated according to the data source identifier into preprocessed multi-source terrain grid data.

[0029] It should be noted that the Hamper filtering rule achieves outlier detection and replacement based on the ranking relationship between the median and absolute deviation of the neighborhood: a neighborhood window is constructed with a certain grid cell in the unified reference grid as the center. The side length of the neighborhood window is determined by rounding up and taking an odd number from the "equivalent number of units in the unified reference grid converted from the resolution field" in the data registration record; the median value of the elevation value in the neighborhood window is taken as the reference value; the absolute deviation between the elevation of each grid cell in the neighborhood window and the reference value is calculated, and the absolute deviation of the central grid cell is compared with the absolute deviation of the other grid cells in the neighborhood window one by one; when the absolute deviation of the central grid cell reaches the maximum value in the neighborhood window, the central grid cell is marked as an outlier, and the elevation of the central grid cell is replaced with the reference value, thereby achieving outlier detection and replacement without introducing additional thresholds or coefficients.

[0030] S2: Based on the preprocessed multi-source terrain grid data, the slope, curvature and landform boundary line are calculated to form a structure guidance field. A confidence model is established based on sensor error, spatial density and landform consistency score to generate a confidence field. A structure-preserving fusion algorithm is used to generate a unified terrain field. S2.1: Extract source elevation grids, overlay marker grids, missing measurement masks, and filling masks from the preprocessed multi-source terrain grid data, and establish an effective cell set; Furthermore, the source elevation grid, overlay marker grid, missing measurement mask, and filling mask are read according to the data source identifier and aligned according to the unified reference grid index. Within each data source, the grid cell index with an overlay marker value of 1 and a missing measurement mask value of 0 is written into the valid cell set, and the grid cell index with a filling mask value of 1 is written into the filling cell set. At the multi-source level, the valid cell set is counted according to the grid cell index, and the grid cell indexes of data sources with a count of more than two are written into the overlapping overlay grid cell set.

[0031] S2.2: Calculate the slope and curvature within the effective cell set, generate slope and curvature rasters, construct boundary response rasters on a unified reference grid, and extract the terrain boundary line mask; Furthermore, within the overlapping grid cell set, the median elevation of the multi-source elevations of the same grid cell is taken and written into the reference elevation raster; the slope raster is calculated using central difference on the reference elevation raster, and the curvature raster is calculated using Laplace second-order difference; when constructing the boundary response raster, the range normalization of the neighborhood difference amplitude of the slope raster, the neighborhood difference amplitude of the curvature raster, and the dispersion of the multi-source elevations of the same grid cell are respectively summed and written back; the boundary response raster is segmented using Otsu's method and written into the candidate boundary zone; the candidate boundary zone is filtered by connected component to remove isolated patches, and then morphological refinement is used to write into the landform boundary line mask.

[0032] The slope grid is calculated using the central difference method, and the expression is: ; in, In the grid cell The slope grid at that location, In the grid cell Reference elevation grid to the east, In the grid cell Reference elevation grid to the west, In the grid cell Reference elevation grid to the north, In the grid cell Reference elevation grid in the south, To standardize the east-west grid step size of the reference grid, To standardize the north-south grid step size of the reference grid; It should be noted that the reference elevation raster is a numerical raster that uses a unified reference grid as its carrier and assigns a "baseline elevation value" to each grid cell under the same grid indexing system. It is used as a comparison benchmark and input for morphological calculations after multi-source elevation alignment. Its generation process involves taking a robust representative value of the multi-source elevation of the same grid cell within the overlapping grid cell set and writing it into that grid cell. This ensures that the reference elevation raster is consistent with subsequent processes such as slope raster, curvature raster, error benchmark, and offset identification in terms of spatial range, grid step size, and grid index. This supports operations such as residual calculation, boundary response construction, and structural constraint fusion iteration for various source elevation grids.

[0033] The curvature raster is calculated using the second-order Laplace difference, and the expression is: ; in, In the grid cell Curvature grid at the location, In the grid cell Reference elevation grid at the location; It should be noted that Otsu's threshold segmentation method is an automatic threshold selection method based on histograms. First, the values ​​of the raster to be segmented are statistically converted into a histogram. All possible thresholds are used as segmentation points to divide the pixels into two categories: "below the threshold" and "above the threshold". For each candidate threshold, the pixel proportion and intra-class mean of the two categories are calculated, and the inter-class variance is calculated accordingly. The threshold that maximizes the inter-class variance is selected as the segmentation threshold. The raster is converted into a binary result using the segmentation threshold, thereby achieving automatic segmentation of the foreground and background without relying on manual threshold setting.

[0034] It should be noted that Otsu's threshold segmentation method is an automatic threshold selection method based on histograms. First, the values ​​of the raster to be segmented are statistically converted into a histogram. All possible thresholds are used as segmentation points to divide the pixels into two categories: "below the threshold" and "above the threshold". For each candidate threshold, the pixel proportion and intra-class mean of the two categories are calculated, and the inter-class variance is calculated accordingly. The threshold that maximizes the inter-class variance is selected as the segmentation threshold. The raster is converted into a binary result using the segmentation threshold, thereby achieving automatic segmentation of the foreground and background without relying on manual threshold setting.

[0035] S2.3: Combine the slope grid, curvature grid, and terrain boundary mask into a structure guiding field according to a unified reference grid index; Furthermore, verify that the slope raster, curvature raster, and landform boundary mask are consistent in spatial range, grid step size, and grid index system; write the slope raster to the slope layer, the curvature raster to the curvature layer, and the landform boundary mask to the boundary constraint layer; write the three layers side by side to the same grid cell position in the structure guidance field according to the unified reference grid index.

[0036] It should be noted that the slope layer uses the slope raster value of each grid cell as the attribute of that cell, and is used to characterize the strength of local undulations and the sensitive area of ​​slope direction changes; the curvature layer uses the curvature raster value of each grid cell as the attribute of that cell, and is used to characterize the topographic undulation and the strength of bends; the boundary constraint layer uses the binary label of the landform boundary line mask as the attribute of that cell, and marks whether the cell is located at the landform boundary line, so as to apply isolation constraints to cross-boundary diffusion, interpolation neighborhood selection or update propagation in subsequent structure-preserving fusion, resampling and smoothing updates; the three are stored side by side at the same grid cell location to form a structure guiding field, so that subsequent calculations can read the slope, curvature and boundary label at the same index and work together to constrain the data.

[0037] S2.4: Construct error baseline quantities on a unified reference grid and convert them into error confidence through a monotonically decreasing mapping to generate an error confidence layer; Furthermore, for each data source, the residual magnitude of the source elevation grid relative to the reference elevation raster is calculated within the effective cell set. The median value of the residual magnitude is then written into the residual statistics layer using a neighborhood window. The side length of the neighborhood window is rounded up to an odd number based on the number of equivalent cells on the unified reference grid in the resolution field. The nominal error level in the data registration record is mapped to the nominal error amount. The residual statistics layer and the nominal error layer are added together and written into the error benchmark layer. The minimum and maximum values ​​within the effective cell set of the error benchmark layer are linearly scaled to obtain the normalized error amount. The normalized error amount is then reverse-scaled and written into the error confidence layer. Invalid cells are marked as invalid.

[0038] It should be noted that the error benchmark layer is a raster layer that uses a unified reference grid as a carrier to quantify and record the "elevation error level of each grid cell". The grid cell values ​​are derived from the merging of two parts of information: one part comes from the residual statistics obtained by the data source relative to the reference elevation raster within the effective cell set (used to reflect the actual deviation level at that location), and the other part comes from the nominal error of the data source corresponding to the data source in the data registration record (used to reflect the inherent measurement accuracy of the data source). After the two parts are merged and written under the same grid index, the error benchmark layer spatially possesses two types of constraint information: "local observation deviation" and "data source accuracy level". When it is subsequently converted into an error confidence layer through a monotonically decreasing mapping, the confidence level of the grid cell with the larger error is lower, thus participating in the construction and fusion of confidence field constraints.

[0039] The residual magnitude of the source elevation grid relative to the reference elevation raster is calculated using the following expression: ; in, For the dimensionless quantity of the residual amplitude, The data source identifier for the source elevation grid is... Time grid cell The elevation value at that location, For reference elevation raster in grid cells The elevation value at that location, The nominal error is the result of mapping the nominal error level to the nominal error in the data registration record.

[0040] S2.5: Calculate the density baseline quantity on the unified reference grid and convert it into density confidence through a monotonically increasing mapping to generate a density confidence layer; Furthermore, for each data source, the local coverage ratio is calculated using a neighborhood window on a unified reference grid. The local coverage ratio is the ratio of the number of cells with a coverage marker value of 1 in the neighborhood to the total number of neighboring cells, and is written into the density reference layer. The side length of the neighborhood window is rounded up to an odd number based on the equivalent number of cells in the resolution field on the unified reference grid. For grid cells with a missing mask value of 1, an invalid marker is written into the density reference layer. For grid cells with a filling mask value of 1, the deviation amplitude between the filling elevation and the median value of the observed elevation in the neighborhood is calculated. The minimum and maximum values ​​of the deviation amplitude are linearly scaled within the filling cell set, and the scaling result is then reverse-scaled and written into the filling consistency value. The product of the local coverage ratio and the filling consistency value is written back into the density reference value of that grid cell. The minimum and maximum values ​​of the valid cells in the density reference layer are linearly scaled and written into the density confidence layer. Invalid cells retain the invalid marker.

[0041] The local coverage ratio is calculated using a neighborhood window on a unified reference grid, expressed as follows: ; in, The data source identifier for the density reference layer is... Time grid cell The local coverage ratio at that location For grid cells The number of grid cells with a cover marker value of 1 within the neighborhood window centered on the target. This represents the total number of grid cells within the neighborhood window.

[0042] S2.6: Based on the structure-guided field, calculate the morphological difference quantity on the unified reference grid, calculate the boundary coincidence mark in the grid cell covered by the landform boundary line mask, combine the morphological difference quantity and the boundary coincidence mark into a consistency reference quantity, and obtain the consistency confidence through the difference reduction mapping, and generate a consistency confidence layer. Furthermore, for each data source, the morphological difference is calculated within the effective cell set. The morphological difference is the sum of the difference between the slope raster calculated from the slope raster of the data source and the slope raster calculated from the reference elevation raster, and the difference between the curvature raster calculated from the curvature raster of the data source and the curvature raster calculated from the reference elevation raster. The result is then smoothed back using neighborhood median filtering. Within the grid cell where the terrain boundary mask value is 1, the boundary response raster of the data source is constructed based on the source elevation grid of the data source, and the boundary band of the data source is written using Otsu's threshold segmentation. When the boundary band of the data source is 1 and the terrain boundary mask value is 1 within the same grid cell, the boundary coincidence marker is written as 1. The consistency benchmark is written as the morphological difference, and within the grid cell where the terrain boundary mask value is 1 and the boundary coincidence marker value is 0, the consistency benchmark is replaced with the maximum value of the consistency benchmark within its neighborhood window. The minimum and maximum values ​​within the effective cell set of the consistency benchmark layer are linearly scaled to obtain the normalized difference, and the normalized difference is inversely scaled and written into the consistency confidence layer.

[0043] S2.7: The error confidence layer, density confidence layer, and consistency confidence layer are weighted and summed according to the combination coefficients to generate a comprehensive confidence score, and then normalized to obtain the confidence field. Furthermore, for each data source, the error confidence layer, density confidence layer, and consistency confidence layer are read on a unified reference grid. The combination coefficients are weighted equally (the combination coefficients are non-negative and normalized within the same grid cell to make the contributions of the three confidence levels to the overall confidence level comparable). The three confidence levels are added together and written into the overall confidence level layer. The overall confidence levels of all data sources are summed within the same grid cell, and the overall confidence level of each data source is divided by this sum and written into the confidence level field. Grid cells with invalid overall confidence levels are not included in the summation of that grid cell.

[0044] To verify that the confidence field can characterize the reliability of the data source elevation at the grid cell level, the "maximum weight of the confidence field" is used as the representative confidence level of the grid cell. This representative confidence level is then subjected to binning statistics, and within the same bin, the bin-average absolute error and the bin-P90 absolute error are calculated separately, thus establishing a correspondence between "bin-average confidence level and error statistic." This correspondence is used to characterize the degree to which the confidence level exhibits a monotonic trend with increasing error and the controllability of the high error tail, such as... Figure 2 As shown.

[0045] S2.8: Perform weighted synthesis on the source elevation grid based on the confidence field to obtain an initial unified topographic field, and use a structure-preserving fusion algorithm to perform fusion iteration to generate a unified topographic field; Furthermore, the elevation grids and confidence fields of each data source are read on the unified reference grid. The elevations of each data source within the same grid cell are weighted and summed according to their confidence weights and written into the initial unified topographic field. Anisotropic diffusion iterative updates are performed on the initial unified topographic field under the constraints of the structure-guided field. Grid cells with a boundary constraint layer value of 1 suppress cross-boundary diffusion, while non-boundary grid cells adjust the diffusion intensity according to the slope layer and curvature layer. In each iteration, the maximum absolute difference between the unified topographic field before and after the update is calculated. The maximum absolute difference is compared with the elevation resolution corresponding to the resolution field. The iteration ends and the unified topographic field is written back when the maximum absolute difference does not exceed the elevation resolution.

[0046] It should be noted that the structure-preserving fusion algorithm refers to a "controlled iterative smoothing-correction" process performed on an initial unified topographic field. Its core operation is not simple global smoothing, but rather restricting the update amount of each grid cell to a range of direction and intensity consistent with the topographic structure: Local gradient and curvature-related changes in the initial unified topographic field are calculated on a unified reference grid, and an adjustment factor for diffusion intensity is constructed for each grid cell; in each iteration, elevation differences are collected from adjacent cells for each grid cell, scaled according to the adjustment factor, accumulated as the update amount, and written back, thus allowing moderate diffusion in flat areas to suppress noise, while diffusion in structurally significant areas such as ridges, valleys, and cliffs is suppressed to avoid being "flattened"; during the iteration process, the maximum absolute difference between the unified topographic field before and after the update is recorded, and the comparison of this monitored quantity with the elevation resolution triggers a stop, ensuring that the number of iterations does not depend on subjective settings and that the update amplitude falls within a resolvable range.

[0047] The structural guidance field constraint is a collective expression of the "spatial rules and boundary rules" of the above-mentioned iterative update. The structural guidance field contains three types of information: slope layer, curvature layer, and boundary constraint layer. The slope layer and curvature layer are used to provide the identification criteria for "whether the location belongs to the terrain structure sensitive area" at the grid cell level, and adjust the diffusion intensity and update direction preference accordingly, so that the update tends to spread along the contour line direction and suppresses cross-slope propagation. The boundary constraint layer is used to provide "prohibited terrain boundaries". When the boundary constraint layer is set to 1 in a certain grid cell, the contribution of the elevation difference of the adjacent pairs across the boundary is directly suppressed or set to zero during the update, so that the updates of the grid cells on both sides do not affect each other, thereby avoiding boundary drift and structural crosstalk. The above constraints work together to ensure that the fusion iteration reduces noise and fills in local inconsistencies while maintaining the continuity of the terrain boundary position and local morphology.

[0048] To demonstrate that the structure-preserving fusion algorithm has an executable iterative termination condition and can converge stably under the constraints of the structure-guided field, the maximum absolute difference of the unified topographic field before and after each iteration is recorded as a convergence monitoring metric. The convergence monitoring metrics of the structure-preserving fusion iteration, the local reconstruction iteration, and the overall smoothing iteration are compared and displayed. This convergence curve characterizes the process of the iteration update amplitude decreasing with each iteration and stopping the iteration when the "maximum absolute difference does not exceed the elevation resolution," as shown below. Figure 3 As shown.

[0049] S3: Based on the data differences between the unified topographic field and the preprocessed multi-source topographic grid data, and combined with the topographic geometric continuity, identify local offset areas with slope anomalies, gradient abrupt changes, and boundary mismatches. Calculate the offset error field using local registration and gradient consistency constraints, and perform error compensation and local reconstruction on the unified topographic field to obtain the corrected topographic field. S3.1: Calculate the difference grids between the source elevation grid and the uniform topographic field within the grid cells marked as valid by the overlay marker grid, and aggregate them into a composite difference grid; Furthermore, under the unified reference grid index system, the source elevation grid and overlay marker raster of each data source are read, and the residual magnitude between the source elevation of the data source and the corresponding unit elevation of the unified topographic field is calculated in the grid cell with an overlay marker value of 1. The residual magnitude is written into the difference raster of the data source. Robust aggregation is performed on the difference raster of multiple data sources according to the grid cell index within the set of overlapping overlay grid cells. The robust aggregation uses the median value. The aggregation result is written into the comprehensive difference raster, and the dispersion of the multi-source residuals is written into the dispersion marker layer.

[0050] S3.2: Calculate the baseline slope raster and baseline gradient raster based on the unified topographic field, and combine them with the comprehensive difference raster to identify slope anomaly indicators, gradient abrupt change indicators and boundary mismatch indicators, and generate a set of offset risk layers; Furthermore, a reference slope raster is calculated using central difference on a unified topographic field, and a reference gradient raster is calculated using discrete gradient magnitude on the same topographic field. The intensity of abrupt differences is calculated using a neighborhood window on the comprehensive difference raster and written into a difference abruptness layer. The reference slope raster and the difference abruptness layer are jointly judged under the same grid index, and grid cells that meet the criteria of "slope in a high-value area and difference abruptness intensity in a high-value area" are written into a slope anomaly indicator layer. The neighborhood difference magnitude of the reference gradient raster is written into a gradient abruptness indicator layer. Within grid cells where the landform boundary mask value is 1, the consistency of the sign of the comprehensive differences on both sides of the boundary and the consistency of the gradient direction are compared, and grid cells that do not meet the consistency criteria are written into a boundary mismatch indicator layer. The slope anomaly indicator, gradient abruptness indicator, and boundary mismatch indicator are written into the offset risk layer set in parallel.

[0051] The formula for calculating the baseline gradient grid is: ; in, As a reference gradient raster in the grid cell The gradient magnitude metric at that point.

[0052] S3.3: Perform thresholding and connected component aggregation on the set of offset risk layers to generate a local offset region mask; Furthermore, within the statistical range defined by the effective cell set, slope anomaly indicator layers, gradient abrupt change indicator layers, and boundary mismatch indicator layers are read respectively. For each indicator layer, its median value and median absolute deviation are calculated, and the threshold for that indicator layer is determined by the sum of the products of the median value, the outlier discrimination coefficient, and the median absolute deviation. The three thresholds are named slope anomaly threshold, gradient abrupt change threshold, and boundary mismatch threshold, respectively. The outlier discrimination coefficient is read from the numerical model grid configuration parameters. If the configuration file is not provided, the outlier discrimination coefficient is written as the coefficient corresponding to three times the median absolute deviation, and a source mark is written in the data registration record. The value range of each threshold is defined as "between the minimum and maximum values ​​of the indicator layer within the effective cell set." When the calculated threshold exceeds... When the specified range is defined, the threshold is clipped to the corresponding boundary and a clipping mark is written. Grid cells whose indicator layer values ​​exceed the corresponding threshold are written to the corresponding binary risk raster. The union of the three types of binary risk raster is taken to obtain a comprehensive risk binary raster, and connected component marking is performed according to eight neighborhoods. Connected component filtering uses a lower limit for area and a lower limit for shape compactness. The lower limit for area and shape compactness are read from the numerical model grid configuration parameters. If the configuration file is not provided, the lower limit for area is written as the area of ​​no less than one unified reference grid cell and no more than the total area of ​​the statistical range. The lower limit for shape compactness is written as a value between zero and one and a source mark is written in the data registration record. Connected components that meet both the lower limit for area and the lower limit for shape compactness are written to the local offset region mask, and those that do not meet the requirements are written with a rejection mark.

[0053] S3.4: Perform local registration within the area covered by the local offset region mask to obtain a set of local displacement vectors, and interpolate the set of local displacement vectors to generate the initial displacement field; Furthermore, within the mask-covered area of ​​the local offset region, control point grids are selected at fixed step sizes, and at each control point, a local window of the same size is extracted from the unified topographic field and the corresponding source elevation grid. For each pair of local windows, the phase correlation method is performed to estimate the planar displacement, and the translation amount corresponding to the position of the correlation peak is written into the local displacement vector set, and the registration confidence is written with the sharpness of the correlation peak. Control points that do not meet the registration confidence requirements are removed. The local displacement vectors of the retained control points are written into the control point displacement table according to the control point coordinates, and inverse distance weighted interpolation is used to interpolate the displacement vector components separately within the mask-covered area. The interpolation results are written into the initial displacement field according to the grid cell index.

[0054] S3.5: The displacement field is updated by applying gradient consistency constraints to generate the offset error field; Furthermore, within the mask-covered area of ​​the local offset region, the initial displacement field is applied to the source elevation grids of each data source for reverse translation sampling, and the displacement-corrected source elevation is written into it. Gradient grids are calculated on the unified topographic field and the displacement-corrected source elevations, and the gradient difference magnitude is written into the gradient inconsistency layer. The displacement field is iteratively updated within the mask-covered area, with the gradient inconsistency layer used as the driving force for displacement correction, and spatial continuity constraints are applied to the displacement field using Laplace smoothing. Zero-flux boundary conditions are applied at the mask boundaries. The maximum change magnitude of the gradient inconsistency layer in two adjacent iterations is used as the convergence requirement as the iteration termination condition, and the converged displacement field is written into the offset error field.

[0055] It should be noted that spatial continuity constraint refers to the smooth update of each grid cell in the displacement field within the mask coverage area of ​​the local offset region, by referring to the displacement values ​​of its neighboring grid cells. This ensures that the displacement changes of adjacent grid cells remain slow and avoids the "bend" phenomenon of sudden increases or decreases in the displacement of a single grid cell. In specific implementation, Laplace smoothing is used to compare the displacement value of each grid cell with the average displacement of its neighborhood, and to slightly pull back the displacement of the grid cell in the direction of the difference, so that the displacement field is spatially continuous.

[0056] Zero-flux boundary conditions refer to restricting the smooth update of the displacement field at the mask boundary in the local offset region to prevent it from propagating across the mask boundary. This ensures that the displacement correction inside the mask is not affected by the displacement values ​​outside the mask, and the displacement correction outside the mask is not "driven" by the displacement correction inside the mask. In specific implementation, when performing neighborhood update calculations on the mesh cells at the mask boundary, only the adjacent mesh cells located inside the mask are selected to participate in the neighborhood averaging and pullback calculations. The adjacent mesh cells outside the mask are not included in the update, thus creating a processing effect at the boundary that is "smooth only inside the boundary and does not spread across the boundary".

[0057] S3.6: Perform error compensation on the unified topographic field based on the migration error field, and perform local reconstruction in the local migration area to generate the corrected topographic field; Furthermore, within the masked area of ​​the local offset region, the offset error field is applied to the grid cell coordinates of the unified topographic field. Reverse displacement sampling is performed on the unified topographic field and written into the offset-aligned topography, with the elevation difference before and after alignment written into the error compensation amount. The error compensation amount is smoothed using the neighborhood median, and boundary constraints are applied in the masked area of ​​the landform boundary line to suppress cross-boundary diffusion. Local reconstruction is performed on the compensated topography within the masked area. The local reconstruction uses anisotropic diffusion update under the constraints of the structure-guided field, and the diffusion on both sides of the boundary is isolated by the boundary constraint layer. The local area and the non-masked area are directly written back using the same grid index to obtain the corrected topographic field.

[0058] To illustrate the improvement effect of error compensation and local reconstruction based on the offset error field on the global error distribution, the absolute errors of the initial unified topographic field of control group A, the unified topographic field of control group C, and the corrected topographic field of the experimental group relative to the true topographic field were calculated on a unified reference grid using the ground truth topographic field elevation as a reference. The cumulative distribution curve of the absolute error of the global grid cells was then statistically analyzed. This cumulative distribution was used to compare the cumulative proportions achieved by each group under the same error threshold, as well as the differences in the low error range, thereby characterizing the comprehensive improvement effect of error compensation and local reconstruction on the error tail and main range. Figure 3 As shown.

[0059] S4: Based on the numerical model grid structure, perform model grid adaptation and terrain resampling based on structure-guided field constraints on the corrected terrain field to generate a grid terrain field. Then, after performing slope restriction strategy on steep areas in the grid terrain field, obtain the overall smoothed terrain field through overall smoothing processing. S4.1: Collect numerical model grid configuration parameters, establish numerical model grid structure, and generate model grid index set; Furthermore, the numerical model grid configuration file and runtime parameter table are read, and the grid type identifier, horizontal projection parameters, vertical datum, horizontal resolution parameters, grid origin, grid direction parameters, boundary clipping rules, and nesting level parameters are extracted and written into the grid configuration record. The grid structure is established according to the grid type identifier. For regular grids, the grid cell center coordinate sequence is written by expanding the grid cells by rows and columns. For hexagonal grids, the center coordinates and adjacency relationships are written according to the cellular topology. For nested grids, the center coordinates of each layer of grid cells are written according to the hierarchical relationship, and the parent-child mapping is written. The grid cell identifier, center coordinates, and adjacency index are written into the model grid index set. The consistency of coordinate units and the integrity of boundary coverage are checked and checked by writing a check mark. At the same time, the slope threshold parameter and the smoothing termination threshold parameter are extracted from the grid configuration record and written into the parameter field.

[0060] It should be noted that the slope threshold parameter is primarily obtained from the "numerical model grid configuration parameters": Parameters related to terrain slope limitations are retrieved from the grid configuration file and runtime parameter table, and written into the "slope threshold parameter" field in angle form. In scenarios where this parameter is not provided in the configuration file, the slope raster is calculated based on the corrected terrain field and grid adjacency relationships. The "maximum allowable slope angle" is used as the threshold, set to 30° as the default recommended value and written to the grid configuration record. Fine-tuning is allowed within the range of 20° to 40° based on the domain's slope statistical distribution to balance numerical stability and geomorphic feature preservation. The meaning of this 30° recommended value and threshold is given in the user manual and source code of the local terrain smoothing script.

[0061] The smoothing termination threshold parameter is used as the iteration stopping criterion, and its value is obtained based on the "model elevation resolution": The elevation resolution quantization step size is read from the numerical model grid configuration parameters and recorded as the "elevation resolution". The "smoothing termination threshold parameter" is set as a proportional threshold of the "elevation resolution" and written into the parameter field. The proportional threshold is set to half the elevation resolution to the elevation resolution, ensuring that iteration stops when the "convergence monitoring quantity" decreases to the point where it no longer causes discernible changes in the topographic field, thus avoiding over-smoothing that introduces geomorphic structure loss. If the grid configuration parameters do not explicitly provide the elevation resolution quantization step size, the coarser of the elevation unit quantization precision in the data registration record and the effective precision of the topographic variables output by the numerical model is used as the "elevation resolution" and written into the record before being assigned a value.

[0062] S4.2: Based on the corrected terrain field, perform terrain resampling on the pattern grid index set to obtain the grid terrain field, and use the structure guiding field as a constraint to obtain the grid terrain field that fits the pattern grid. Furthermore, at the center coordinates of each grid cell given by the pattern grid index set, spatial interpolation sampling is performed from the corrected topographic field and written into the grid topographic field. Regular grids use bilinear interpolation, hexagonal grids use natural neighborhood interpolation, and nested grids use interpolation consistent with their grid type in the sub-grids and consistent resampling at the parent-child boundary. The boundary constraint layer and slope layer in the structure-guided field are read, and same-side constraints are applied to the interpolation neighborhood at the projection position covered by the landform boundary mask, so that the sample points participating in the interpolation are located on the same side of the boundary. Anisotropic weights are applied to the interpolation neighborhood in the high-value region of the slope layer, so that the weights are concentrated along the contour line direction and cross-slope mixing is suppressed. The constraint sampling results are written back to the grid topographic field, and a sampling quality mark is written in each grid cell.

[0063] S4.3: Calculate the slope grid based on the grid terrain field, generate a steep region mask, and perform connected component aggregation to obtain a set of steep regions; Furthermore, the slope raster is calculated on the grid topographic field according to the grid adjacency relationship. The regular grid adopts the center difference, and the hexagonal grid constructs the discrete gradient magnitude by the elevation difference and center distance of the adjacent cells. The nested grid is calculated in each layer according to the adjacency relationship within the layer and the adjacency is filled in at the layer boundary according to the parent-child mapping. The grid cells of the slope raster that are greater than the slope threshold are written into the steep region mask. The steep region mask is marked with connected components according to the grid adjacency relationship. The set of grid cell identifiers of each connected component is written into the steep region set, and connected components with insufficient area are removed.

[0064] The slope raster is calculated based on the grid topographic field, using the following expression: ; in, For the grid terrain field in adjacent grid cell pairs ( The slope grid on the ) For grid topographic field in grid cell The elevation value at that location, For grid topographic field in grid cell The elevation value at that location, For grid cells in the pattern grid index set and Distance between center coordinates.

[0065] It should be noted that the minimum distance between the centers of adjacent grid cells in the pattern grid index set is used as the spatial scale benchmark, and the elevation resolution quantization step size in the grid configuration record is used as the elevation scale benchmark. The "allowable adjacent elevation variation" is set to 1 to 3 times the elevation resolution quantization step size. Then, the spatial scale and elevation scale are converted into angular slope limits as slope thresholds, with a value range of 15 degrees to 40 degrees.

[0066] S4.4: Execute a slope constraint strategy within a set of steep regions to generate a slope-constrained grid terrain field; Furthermore, within each steep region's connected domain, the boundary grid cells of that connected domain are used as constraint boundaries, and the process is traversed inward along the adjacency relationship. For any pair of adjacent grid cells, the center distance is calculated, and the allowable elevation difference is calculated based on the slope threshold parameter. Pairs of adjacent grid cells whose actual elevation difference exceeds the allowable elevation difference are clipped and written back. The clipping method is to keep the average elevation of the two cells unchanged and reduce the difference to the allowable range. When the same grid cell involves clipping in multiple directions, the median value of the candidate elevations obtained in each direction is written back. The connected domain is traversed cyclically until all adjacent grid cell pairs in the connected domain satisfy the allowable elevation difference constraint. The updated grid topography field is written back as the slope-constrained grid topography field, and a constraint mark is written to the grid cells that have been clipped.

[0067] It should be noted that elevation difference constraint refers to setting an "upper limit" for the elevation difference between any two adjacent grid cells within the connected domain corresponding to a set of steep regions, ensuring that the elevation difference of the adjacent pair does not exceed the allowable range determined by the slope threshold parameter and the center distance between the two cells. In practice, the slope threshold parameter is first used as the "maximum allowable slope" criterion, and then the maximum allowable elevation difference of the adjacent cell pair is calculated in combination with the center distance between the adjacent grid cells. The actual elevation difference is compared with the maximum elevation difference. If it exceeds the limit, the elevation difference of the adjacent cell pair is clipped and written back, so that the clipped elevation difference returns to the allowable range. This constrains the elevation changes of adjacent cells pair by pair within the local connected domain, avoiding the formation of excessively steep gradients and reducing the risk of gradient instability in the initial field of the numerical model.

[0068] S4.5: Based on the slope-constrained grid topographic field, the global grid is smoothly updated with the goal of stabilizing the initial gradient of the model, generating an overall smoothed topographic field; Furthermore, a smoothing iteration is performed on the slope-constrained grid topography field. Smoothing adopts anisotropic diffusion and uses the boundary constraint layer of the structure-guided field as an isolation condition. Zero-flux boundary treatment is applied at the location covered by the boundary constraint layer. In each iteration, the gradient raster is calculated in the global grid, and the maximum absolute difference between the gradient raster in this round and the gradient raster in the previous round is calculated and written into the convergence monitoring quantity. The convergence monitoring quantity not exceeding the smoothing termination threshold parameter is used as the iteration termination condition. The iterated grid topography field is written back as the overall smoothed topography field, and an update count mark and gradient stability mark are written to each grid cell.

[0069] In summary, this invention establishes a numerical model terrain processing workflow for geospatial industrial data processing, enabling the fusion and correction of multi-source terrain data under the premise of structural consistency, continuity and stability, and completing model grid adaptation. This results in higher overall consistency of the output terrain field in terms of landform boundaries, morphological details and gradient stability, and lower risk of local offset and mismatch, thereby improving the quality of the initial field of the numerical model and the stability and usability of the simulation results.

[0070] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A numerical model terrain processing method based on multiple fused data, characterized in that: include, Collect multi-source terrain data, preprocess it, and then map it to a unified reference grid to generate preprocessed multi-source terrain grid data; Based on the preprocessed multi-source terrain grid data, slope, curvature and landform boundary line are calculated to form a structure guidance field. A confidence model is established based on sensor error, spatial density and landform consistency score to generate a confidence field. A structure-preserving fusion algorithm is then used to generate a unified terrain field. Based on the data differences between the unified topographic field and the preprocessed multi-source topographic grid data, and combined with the topographic geometric continuity, local offset regions with slope anomalies, gradient abrupt changes, and boundary mismatches are identified. The offset error field is calculated using local registration and gradient consistency constraints, and error compensation and local reconstruction are performed on the unified topographic field to obtain the corrected topographic field. Based on the numerical model grid structure, the corrected topographic field is subjected to model grid adaptation and topographic resampling based on structure-guided field constraints to generate a grid topographic field. After applying a slope constraint strategy to steep areas in the grid topographic field, the overall smoothed topographic field is obtained through overall smoothing processing.

2. The numerical model terrain processing method based on multiple fused data as described in claim 1, characterized in that: The multi-source terrain data includes digital elevation model data, laser point cloud data, remote sensing inversion terrain data, and surveying vector terrain data.

3. The numerical model terrain processing method based on multiple fused data as described in claim 2, characterized in that: The specific steps for generating the preprocessed multi-source terrain grid data are as follows. According to the data registration records, the multi-source terrain data is subjected to unified coordinate benchmark transformation and elevation benchmark conversion, and noise identification and filtering are performed to generate preprocessed multi-source terrain data; The preprocessed multi-source terrain data is mapped to a unified reference grid to generate source elevation grids and overlay marker gratings. Terrain gap filling and outlier detection are performed on the source elevation grid to generate preprocessed multi-source terrain grid data.

4. The numerical model terrain processing method based on multiple fused data as described in claim 3, characterized in that: The specific steps for forming the structural guiding field are as follows. Extract source elevation grids, overlay marker gratings, missing measurement masks, and filling masks from the preprocessed multi-source terrain grid data, and establish an effective cell set; Within the effective cell set, slope and curvature are calculated to generate slope and curvature rasters. Boundary response rasters are constructed on a unified reference grid to extract the terrain boundary line mask. The slope grid, curvature grid, and terrain boundary mask are combined into a structure guidance field using a unified reference grid index.

5. The numerical model terrain processing method based on multiple fused data as described in claim 4, characterized in that: The specific steps for establishing a confidence model based on sensor error, spatial density, and terrain consistency score are as follows. An error baseline quantity is constructed on a unified reference grid and converted into an error confidence level through a monotonically decreasing mapping, thereby generating an error confidence level layer. Density baseline quantities are calculated on a unified reference grid and converted into density confidence values ​​through a monotonically increasing mapping, generating a density confidence layer. Based on the structure-guided field, the morphological difference quantity is calculated on a unified reference grid, and the boundary coincidence mark is calculated on the grid cell covered by the landform boundary line mask. The morphological difference quantity and the boundary coincidence mark are combined into a consistency benchmark quantity and the consistency confidence is obtained through difference reduction mapping, and a consistency confidence layer is generated.

6. The numerical model terrain processing method based on multiple fused data as described in claim 5, characterized in that: The specific steps for generating a unified terrain field are as follows. The error confidence layer, density confidence layer, and consistency confidence layer are weighted and summed according to the combination coefficients to generate a comprehensive confidence score, which is then normalized to obtain the confidence field. The source elevation grid is weighted and synthesized based on the confidence field to obtain an initial unified topographic field. A structure-preserving fusion algorithm is then used to perform fusion iterations to generate a unified topographic field.

7. The numerical model terrain processing method based on multiple fused data as described in claim 6, characterized in that: The specific steps for identifying local offset regions with slope anomalies, gradient abrupt changes, and boundary mismatches are as follows. Within the grid cells marked as valid by the overlay marker raster, calculate the difference raster between the source elevation grid and the uniform topographic field, and aggregate them into a comprehensive difference raster; The baseline slope raster and baseline gradient raster are calculated based on the unified topographic field, and combined with the comprehensive difference raster, slope anomaly indicators, gradient abrupt change indicators and boundary mismatch indicators are identified to generate a set of offset risk layers. Thresholding and connected component aggregation are performed on the set of offset risk layers to generate a local offset region mask.

8. The numerical model terrain processing method based on multiple fused data as described in claim 7, characterized in that: The specific steps for performing error compensation and local reconstruction on the unified topographic field to obtain the corrected topographic field are as follows: Local registration is performed within the area covered by the local offset region mask to obtain a set of local displacement vectors, and the set of local displacement vectors is interpolated to generate the initial displacement field; The displacement field is updated by applying gradient consistency constraints to generate the offset error field. Error compensation is performed on the uniform topographic field based on the migration error field, and local reconstruction is performed in the local migration area to generate the corrected topographic field.

9. The numerical model terrain processing method based on multiple fused data as described in claim 8, characterized in that: The specific steps for generating the grid terrain field are as follows: Collect numerical model grid configuration parameters, establish numerical model grid structure, and generate a set of model grid indexes; Based on the corrected terrain field, terrain resampling is performed on the pattern grid index set to obtain the grid terrain field, and the structure guiding field is used as a constraint to obtain the grid terrain field that fits the pattern grid.

10. The numerical model terrain processing method based on multiple fused data as described in claim 9, characterized in that: The specific steps for obtaining the overall smoothed terrain field are as follows. Calculate the slope grid based on the grid terrain field, generate a steep region mask, and perform connected component aggregation to obtain a set of steep regions; Apply a slope constraint strategy within a set of steep regions to generate a slope-constrained grid terrain field; Based on the slope-constrained grid topographic field, the global grid is smoothly updated with the goal of stabilizing the initial gradient of the model, resulting in an overall smoothed topographic field.

Images