A method for generating a seabed DEM based on large-scale chart depth points and high-resolution remote sensing surface

Abstract

This invention provides a method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, belonging to the field of marine surveying technology. The method includes: acquiring discrete depth points from remote sensing images and nautical charts; extracting a relative depth trend surface based on a dual-band logarithmic ratio model; calculating the image gradient field and pixel-level structure tensor, solving for the principal eigenvectors, and constructing an anisotropic diffusion tensor by combining local gradient magnitudes; using the nautical chart depth points to perform regression calibration on the trend surface to generate sparse residuals, and establishing a residual calculation grid; using the anisotropic diffusion tensor as a guiding constraint, performing non-uniform diffusion calculation on the sparse residuals to generate a global continuous residual correction surface; superimposing the correction surface onto the trend surface, and obtaining the final seabed DEM through anomaly removal and edge smoothing. This invention solves the problems of low accuracy and loss of detail in complex shallow sea topography inversion.

Description

Technical Field

[0001] This invention relates to the field of marine surveying technology, and in particular to a method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces. Background Technology

[0002] Currently, the construction of seabed digital elevation models mainly relies on shipborne sonar measurements and airborne lidar measurements. Although these methods offer high accuracy, their high operating costs and long operation cycles make it difficult to achieve rapid updates of data over large areas of the sea. In contrast, while large-scale electronic nautical charts contain officially approved precise depth points, their data distribution is discrete and non-uniform with large point spacing, making it difficult to depict seabed topographic features in detail. While high-resolution multispectral remote sensing imagery can provide extensive and continuous seabed spectral texture information, its direct inversion of depth data often fails to meet the absolute vertical accuracy requirements of engineering projects due to limitations in water turbidity, variations in seabed sediment types, and the complexity of radiative transfer models.

[0003] To balance accuracy and efficiency, current technological trends are increasingly shifting towards heterogeneous fusion of multi-source data. This involves using nautical chart depth points as control points and combining them with spectral information from remote sensing imagery for collaborative inversion. Mainstream research directions include nonlinear regression models based on random forests or neural networks, and spatial nonstationarity modeling based on geographic weighted regression. The core idea is to use discrete, precise depth points to correct continuous remote sensing inversion surfaces, thereby obtaining large-scale, high-precision depth data.

[0004] However, existing multi-source fusion methods still have significant drawbacks. First, conventional interpolation methods only consider Euclidean distance and ignore the anisotropy of seabed topography, resulting in the loss of high-frequency details in the generated digital elevation models at abrupt topographic changes such as sand waves and reef edges due to over-smoothing. Second, most existing point-to-area fusion models are based on statistical regression of global or local windows, assuming that the error is spatially uniform or smoothly varied. However, in reality, the error distribution of remote sensing inversion is highly correlated with seabed sediment type and image texture structure. Existing technologies lack a mechanism that can actively guide the non-uniform propagation and correction of water depth errors using the texture structure features of the image itself, resulting in insufficient fidelity of the synthesized digital elevation models in areas with complex topographic textures. Summary of the Invention

[0005] To overcome the shortcomings of existing technologies, the purpose of this invention is to provide a method for generating seabed DEMs based on large-scale nautical chart depth points and high-resolution remote sensing surfaces. This invention solves the problem that existing technologies neglect to utilize image texture structure features to guide the anisotropic propagation of depth inversion errors, resulting in the loss of details in complex terrain areas.

[0006] To achieve the above objectives, the present invention provides the following solution: A method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces includes: High-resolution remote sensing images and discrete water depth points on nautical charts of the target sea area are acquired. Radiometric and geometric corrections are performed on the high-resolution remote sensing images to obtain target images of the water-bound area. Based on a dual-band logarithmic ratio model, a relative water depth trend surface is extracted from the target images of the water-bound area. Gradient field calculation is performed on the target image of the water-contaminated area to obtain the local gradient magnitude and pixel-level structure tensor matrix, and feature decomposition is performed on the pixel-level structure tensor matrix to obtain the main feature vector. An anisotropic diffusion tensor is constructed based on the principal eigenvectors and local gradient magnitudes; Within the effective mask range, the relative water depth trend surface is regressed and calibrated using the discrete water depth points on the nautical chart to generate a sparse residual point set, and a residual calculation grid is established based on the pixel resolution of the target image of the water-intrusion area. The effective mask range is a connected region with continuous water spectral characteristics extracted from the target image of the water-intrusion area using the normalized water index. An initial error field is constructed on the residual calculation grid using the sparse residual point set, and the anisotropic diffusion tensor is used as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a global continuous residual correction surface. The global continuous residual correction surface is superimposed on the relative water depth trend surface to obtain the initial seabed elevation model; Outlier removal and edge smoothing are performed on the initial seabed elevation model to obtain the final seabed DEM.

[0007] The present invention discloses the following technical effects: This invention provides a method for generating a digital elevation model (DEM) based on large-scale nautical chart depth points and high-resolution remote sensing surfaces. By introducing an anisotropic diffusion mechanism guided by image texture structure tensors, this invention effectively solves the problem of smooth transitions in details caused by neglecting topographic anisotropy in conventional interpolation methods. This invention utilizes the texture gradient features of remote sensing images to construct a differentiated diffusion field, forcing the depth inversion residuals to preferentially propagate non-uniformly along topographic texture directions (such as sand ridges and gullies), and blocking error crossings at abrupt topographic changes, thus breaking the limitation of the assumption of uniform error spatial distribution in traditional point-surface fusion models. This method achieves high-precision control of discrete depth points on large-scale nautical charts and deep fusion of texture details from high-resolution remote sensing images, significantly improving the topographic fidelity and vertical inversion accuracy of digital elevation models (DEMs) generated in complex topographic areas such as sand dunes and reefs. Attached Figure Description

[0008] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. 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.

[0009] Figure 1 This is a flowchart illustrating a method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, provided as an embodiment of the present invention. Detailed Implementation

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

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

[0012] like Figure 1 As shown, this invention provides a method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, including: Step 100: Acquire high-resolution remote sensing images and discrete water depth points of the target sea area, perform radiometric and geometric correction on the high-resolution remote sensing images to obtain the target image of the water-bound area, and extract the relative water depth trend surface from the target image of the water-bound area based on the dual-band logarithmic ratio model. Step 200: Calculate the gradient field of the target image of the water-contaminated area to obtain the local gradient magnitude and pixel-level structure tensor matrix, and perform feature decomposition on the pixel-level structure tensor matrix to obtain the main feature vector. Step 300: Construct an anisotropic diffusion tensor based on the principal eigenvectors and local gradient magnitudes; Step 400: Within the effective mask range, the relative water depth trend surface is regressed and calibrated using the discrete water depth points on the nautical chart to generate a sparse residual point set, and a residual calculation grid is established based on the pixel resolution of the target image of the water-bound area. The effective mask range is a connected region with continuous water spectral characteristics extracted from the target image of the water-bound area using the normalized water index. Step 500: Construct an initial error field on the residual calculation grid using the sparse residual point set, and use the anisotropic diffusion tensor as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a global continuous residual correction surface. Step 600: Superimpose the global continuous residual correction surface onto the relative water depth trend surface to obtain the initial seabed elevation model; Step 700: Perform outlier removal and edge smoothing on the initial seabed elevation model to obtain the final seabed DEM.

[0013] Furthermore, the specific implementation process of step 100 is as follows: This embodiment first acquires high-resolution multispectral remote sensing image data of the target sea area and discrete depth point data of a large-scale electronic nautical chart for the region. To eliminate the influence of atmospheric scattering and sensor errors, this embodiment performs radiometric calibration and atmospheric correction on the original remote sensing image, converting it into a physically meaningful surface reflectance image. Based on this, this embodiment calculates the normalized water index using the green and near-infrared bands and selects the optimal land-water segmentation threshold based on histogram statistical characteristics to construct a binary water mask. By performing spatial logical AND operation on the surface reflectance image using this mask, this embodiment can accurately remove non-water-related pixels such as land, clouds, and man-made structures, thereby obtaining a target image of the water-related area containing only water spectral information.

[0014] For the acquired target image of the water-bound area, this embodiment further separates the blue band reflectance matrix and green band reflectance matrix, which are most sensitive to water depth inversion, from the multispectral data cube. This embodiment uses a dual-band logarithmic ratio model as the basic physical model for water depth inversion. The physical mechanism of this model lies in using the difference in the attenuation coefficients of different bands in the water body to offset the interference caused by changes in seabed type. Specifically, this embodiment performs a natural logarithmic transformation on each effective pixel value in the blue band reflectance matrix and the green band reflectance matrix, and then calculates the ratio of the logarithmically transformed blue band value to the green band value. Through this pixel-level algebraic operation, this embodiment generates an original relative water depth matrix that reflects the relative undulation characteristics of the seabed topography.

[0015] Considering that the original relative depth matrix inevitably contains high-frequency salt-and-pepper noise caused by sensor noise or fine sea surface ripples, this embodiment does not use it directly, but instead performs spatial filtering. This embodiment uses a median filtering algorithm to denoise the original relative depth matrix, setting an appropriate filtering window size and performing sliding statistics within the window, replacing the center pixel value with the median value within the window. This processing step effectively smooths random noise while preserving the gradient characteristics of key terrain features such as sand waves and reef edges, ultimately obtaining a relative depth trend surface with a high signal-to-noise ratio and clear texture, laying a high-quality data foundation for subsequent absolute depth calibration using nautical chart depth points.

[0016] Furthermore, the specific implementation process of step 200 is as follows: This embodiment first performs pixel-level gradient calculations on the target image of the water-bound area to capture the spectral texture variation characteristics of the seabed topography in remote sensing imagery. This embodiment uses a preset gradient convolution kernel to perform sliding convolution operations on the image in both the horizontal and vertical directions, thereby obtaining the horizontal and vertical gradient components of the image in a two-dimensional plane. These two components represent the rate of change of pixel grayscale in the row and column directions, respectively. Based on this, this embodiment performs vector magnitude operations on the horizontal and vertical gradient components at each pixel location, specifically by calculating the square root of the sum of the squares of the two components to obtain the local gradient magnitude. This local gradient magnitude intuitively reflects the intensity information of the image texture edges, providing a quantitative basis for subsequent identification of abrupt terrain changes.

[0017] To accurately describe the local geometric structure and orientation of image texture, this embodiment further constructs a pixel-level structure tensor matrix. This embodiment utilizes the gradient components calculated above to construct a gradient autocorrelation matrix. Each element of this matrix consists of the squares of the horizontal gradient components, the squares of the vertical gradient components, and the products of the horizontal and vertical gradient components. Since single-point gradients are often susceptible to noise interference and lack spatial continuity, this embodiment uses a Gaussian convolution kernel to perform weighted spatial smoothing on the gradient autocorrelation matrix. This step is equivalent to performing a weighted integration of the gradient information within the local neighborhood window of each pixel, thereby generating a pixel-level structure tensor matrix that can robustly represent the dominant direction of local texture.

[0018] For the generated pixel-level structure tensor matrix, this embodiment uses algebraic methods to analyze its inherent geometric features. This embodiment performs eigenvalue decomposition on the structure tensor matrix at each pixel, thereby resolving two non-negative eigenvalues ​​and two orthogonal eigenvectors corresponding to these two eigenvalues. Physically, the magnitude of the eigenvalue reflects the degree of change in local image texture along the direction of the corresponding eigenvector. Therefore, this embodiment compares the absolute magnitudes of the two eigenvalues ​​and determines the eigenvector corresponding to the eigenvalue with the larger absolute value as the principal eigenvector. This principal eigenvector indicates the direction of the fastest change in local image texture, i.e., the normal direction perpendicular to the seabed topographic texture, providing crucial directional guidance for the subsequent construction of the anisotropic diffusion tensor.

[0019] Furthermore, the specific implementation process of step 300 is as follows: This embodiment first derives and constructs a secondary eigenvector to describe the texture extension direction based on the principal eigenvector calculated in the previous steps. Since the principal eigenvector physically represents the gradient direction (i.e., the edge normal direction) where local texture changes are most drastic, this embodiment, based on the principle of vector orthogonality in two-dimensional Euclidean space, calculates a unit vector perpendicular to the principal eigenvector and with normalized magnitude through algebraic operations of coordinate component interchanging and sign inversion, and defines this as the secondary eigenvector. Geometrically, this secondary eigenvector represents the tangential extension direction of topographic textures such as seafloor ridges and gullies, and together with the principal eigenvector, constitutes an orthogonal coordinate basis capable of fully describing the anisotropic geometric features of local terrain.

[0020] To achieve intelligent control of the propagation path of water depth inversion residuals, this embodiment needs to determine the diffusion intensity along the normal and tangential directions of the terrain texture, respectively. This embodiment introduces a preset edge stopping function to perform nonlinear mapping calculations on the local gradient magnitude to obtain the normal diffusion coefficient. The logic design of this edge stopping function aims to simulate the physical barrier effect. When a significant increase in the local gradient magnitude indicates the presence of abrupt terrain changes or edges, the normal diffusion coefficient calculated by the function will approach zero, thereby blocking the propagation of errors across edges; conversely, diffusion is allowed in flat areas. Simultaneously, this embodiment sets a relatively large and constant preset parameter as the tangential diffusion coefficient, aiming to force error information to undergo long-distance smoothing and propagation along the continuous direction of the terrain texture, thereby maintaining the continuity of the terrain structure.

[0021] After establishing the orthogonal direction vectors and their corresponding diffusion control parameters, this embodiment constructs an anisotropic diffusion tensor using matrix synthesis techniques. First, a fundamental matrix is ​​generated by multiplying the principal eigenvectors with their transposes. This fundamental matrix is ​​then weighted using the normal diffusion coefficient to construct the normal projection matrix. Similarly, the tangential projection matrix is ​​constructed by multiplying the secondary eigenvectors with their transposes and using the tangential diffusion coefficient. Finally, matrix addition is performed on the normal and tangential projection matrices to fuse the diffusion constraint information from the two orthogonal directions into a single second-order tensor matrix. This anisotropic diffusion tensor serves as the core guiding operator for solving subsequent partial differential equations, determining the non-uniform flow characteristics during residual correction.

[0022] Specifically, the expression for the edge stopping function is: ; in, The normal diffusion coefficient is the output of the edge stopping function; The magnitude of the local gradient; The gradient scale parameter is used to characterize the boundary scale of edge intensity; This is a grayscale or single-band reflectance intensity map of the target image in the water-contaminated area. The spatial gradient operator for image intensity; Let be the Euclidean norm of the vector.

[0023] Specifically, in this embodiment, to achieve intelligent truncation of the water depth inversion residual propagation process, an edge stopping function based on an inverse square decay mechanism is used to dynamically solve the normal diffusion coefficient. The specific calculation logic of this edge stopping function is as follows: First, the ratio between the local gradient magnitude of the current pixel in the target image of the wading area and the preset gradient scale parameter is calculated. Then, the square of this ratio is calculated and added to a constant to construct the decay denominator. Finally, the reciprocal of the decay denominator is taken as the calculation result, which is the normal diffusion coefficient corresponding to the pixel.

[0024] In this calculation, the normal diffusion coefficient is the output variable of the function, with its value strictly constrained between zero and one. Its function is to quantify and control the ability of the residual correction to penetrate across texture edges. A value closer to one indicates a flatter location, allowing errors to penetrate freely; a value closer to zero indicates a steep edge, blocking error propagation. The local gradient magnitude is the input variable of this function, derived from gradient field calculations performed on the target image of the wading area. It characterizes the degree of change in spectral texture or terrain undulation at the current pixel location. The gradient scale parameter is a preset control constant that defines a soft boundary threshold between "edges" and "flat areas," determining the sensitivity of the diffusion process to gradient changes. In practice, the value of this parameter is usually set according to the overall noise level of the image. For example, for remote sensing images with a gray level range of 0 to 255, this parameter can be selected as an empirical value between 10 and 50. When the local gradient magnitude is much smaller than this parameter, the function output approaches 1, which is an isotropic diffusion to smooth the noise. When the local gradient magnitude is much larger than this parameter, the function output decays rapidly and approaches 0, which is a cessation of diffusion to maintain the terrain features.

[0025] More specifically, this embodiment first receives the pixel-level principal feature vector obtained from the feature decomposition in the preceding step. Since this principal feature vector is essentially a geometric direction vector defined in a two-dimensional image plane, this embodiment performs coordinate analysis on it, decoupling it into two independent scalar values ​​to facilitate subsequent algebraic operations. Specifically, this embodiment extracts the projection values ​​of the principal feature vector onto the image row direction coordinate axis as the principal horizontal component value, and simultaneously extracts the projection values ​​of the principal feature vector onto the image column direction coordinate axis as the principal vertical component value. This step aims to convert the encapsulated vector object into a basic data unit that can be directly transformed numerically, providing accurate data input for subsequent vector derivation based on the principle of geometric orthogonality.

[0026] Subsequently, this embodiment determines the specific component values ​​of the secondary eigenvector based on the orthogonal geometric principles of two-dimensional vector space. Since the secondary eigenvector needs to be strictly perpendicular to the principal eigenvector in its geometric definition to construct a locally orthogonal coordinate system, this embodiment employs an algebraic strategy of component interchanging and inverting to calculate the coordinates of the secondary eigenvector. In the specific calculation logic, this embodiment sets the value of the secondary horizontal component to be determined as the negative of the principal vertical component value obtained from the aforementioned analysis, that is, multiplying the principal vertical component value by -1; simultaneously, this embodiment directly sets the value of the secondary vertical component to be determined as equivalent to the principal horizontal component value obtained from the aforementioned analysis. Through this specific numerical mapping relationship, this embodiment mathematically guarantees that the sum of the products of the two generated vector components is always 0, thus satisfying the orthogonal constraint condition that the vector dot product is 0.

[0027] Finally, in this embodiment, the calculated secondary horizontal and secondary vertical component values ​​are recombined according to a standard two-dimensional vector format. This embodiment constructs a new vector data structure, assigning the secondary horizontal component values ​​to coordinate elements in their horizontal dimension and the secondary vertical component values ​​to coordinate elements in their vertical dimension, thereby generating a complete secondary feature vector. This secondary feature vector physically represents the extension direction of the image texture, i.e., the tangential direction. Together with the principal feature vector representing the texture gradient normal, it constitutes a complete local texture description basis, providing necessary directional guidance data for the subsequent construction of the normal projection matrix and the tangential projection matrix.

[0028] Furthermore, the specific implementation process of step 400 is as follows: This embodiment first initiates a process of traversing and filtering the discrete depth point set on the nautical chart to ensure that the control points involved in the calculation are strictly located within areas with water spectral characteristics. This embodiment reads the geographic coordinates of each discrete depth point and projects them onto the raster coordinate system of the target image of the water-bound area, determining whether the coordinate point falls within the effective mask range determined by the normalized water index. For the selected effective depth control points falling within this mask range, this embodiment extracts the absolute depth recorded in its nautical chart attributes as the true depth value, and simultaneously extracts the relative light intensity ratio as the trend depth value based on the corresponding pixel position of the relative depth trend surface according to the planar coordinate index of the point. Through this precise matching of spatial locations, this embodiment establishes a series of regression sample pairs containing dimensionless relative depth and true physical depth.

[0029] Based on the constructed regression sample pair sequence, this embodiment employs the least squares method for linear regression analysis to solve for the water depth inversion parameters. This embodiment calculates the optimal conversion gain coefficient and conversion offset coefficient by minimizing the sum of squares of the differences between the actual and predicted water depths at all sample points, thereby establishing a linear water depth inversion model from the relative trend surface to the absolute water depth surface. Using this model, this embodiment converts the previously extracted trend depth values ​​into theoretical model-predicted water depth values, and subtracts the model-predicted water depth value from the corresponding actual water depth value to calculate the inversion residual value for each control point location. This inversion residual value physically quantifies the nonlinear terrain details or local biases that the linear model fails to explain. This embodiment combines and stores the coordinates of all effective control points with their corresponding inversion residual values, forming a sparse residual point set.

[0030] To provide a standard numerical space carrier for subsequent residual diffusion calculations, this embodiment constructs a residual calculation grid strictly aligned with the image space. This embodiment directly reads the metadata information of the target image of the water-affected area, obtaining its matrix row count, matrix column count, and geospatial resolution parameters. Based on these spatial dimension parameters, this embodiment initializes a two-dimensional matrix in memory with the exact same size as the target image of the water-affected area, setting the initial value of all elements in this matrix to zero. This embodiment marks this zero-value matrix as the residual calculation grid, where each grid node corresponds to a resolution cell in the actual sea area, thus providing a basic computational domain for extending sparsely distributed residual data to the entire domain.

[0031] Furthermore, the expression for the linear water depth inversion model is: ; in, coordinates The model predicts the water depth at that location; For the relative water depth trend surface in Trend depth value at the location; The conversion gain coefficient is used to map the trend depth scale to the water depth scale; The conversion offset coefficient is used for reference translation; These are the geographic spatial plane coordinates.

[0032] Specifically, this embodiment employs a linear regression mathematical model to establish the numerical mapping relationship between dimensionless relative spectral depth and physically meaningful absolute water depth. The specific calculation logic of this model is configured as follows: the trend depth value extracted from the relative water depth trend surface at the current geospatial plane coordinates is used as the input independent variable. First, it is multiplied by a defined conversion gain coefficient. Then, the result of this multiplication is added to a defined conversion offset coefficient, and finally, the model-predicted water depth value at that coordinate is output.

[0033] In this calculation process, the model-predicted water depth value is the output target of the algorithm, referring to the vertical water depth data in meters after correction by regression parameters; the trend depth value is the input basis of the algorithm, referring to the relative topographic relief value calculated using the logarithmic ratio of the blue-green bands of remote sensing images. The conversion gain coefficient is a scaling factor used for scale adjustment. It is derived from the least squares linear fit of the nautical chart control point pairs and functions to map the dimensionless relative rate of change of the spectrum to the physical rate of change of the actual water depth. For example, in practical implementation, if regression analysis shows that each unit increase in relative spectral value corresponds to an increase of 25.5 meters in water depth, then the coefficient is determined to be 25.5; the transformation offset coefficient is an intercept constant used for benchmark correction, also derived from least squares fitting, and its function is to eliminate systematic fixed errors caused by tidal height differences or atmospheric correction residuals, and to perform vertical translation of the reference plane. For example, if the intercept of the fitted line is -1.2, then the coefficient is set to -1.2; the geospatial plane coordinates are used to uniquely identify the currently calculated pixel position on the projected map.

[0034] Furthermore, the specific implementation process of step 500 is as follows: This embodiment first performs a discrete mapping operation from sparse data to a continuous field space to construct an initial state field for numerical computation. This embodiment reads the pixel coordinate index and corresponding inversion residual value of each point in the sparse residual point set and accurately assigns them to the corresponding pixel positions in the residual computation grid. For blank pixel positions in the grid that do not have direct assignments, this embodiment initializes their values ​​to zero, thus forming an initial error field. Simultaneously, to ensure that the accuracy of known control points is not compromised by the smoothing effect during subsequent diffusion, this embodiment records and marks all pixel positions mapped to true inversion residual values ​​as a fixed source term coordinate set. This coordinate set will serve as the Dirichlet boundary condition or mandatory constraint term in subsequent iterative calculations.

[0035] Subsequently, this embodiment enters the numerical iterative solution stage for non-uniform diffusion, the core of which lies in using a tensor field to guide the direction of error flow. In each iteration cycle, this embodiment first uses a difference operator to perform gradient calculation on the current initial error field, obtaining an error gradient field characterizing the rate of change of the error space. Next, this embodiment performs pixel-by-pixel matrix multiplication between the anisotropic diffusion tensor constructed in the previous steps and the error gradient field, generating a diffusion flux vector field. In this process, the anisotropic diffusion tensor acts as a conduit, forcing the error flux to flow preferentially along the tangential direction of the image texture and suppressing it in the texture normal, i.e., the edge direction. Subsequently, this embodiment performs divergence calculation on the diffusion flux vector field to obtain the error update matrix, and combines this with a preset time step to numerically accumulate and update the initial error field, thereby obtaining the iterative intermediate field at that time step.

[0036] After updating the time step in a single iteration, this embodiment executes forced constraints and convergence determination logic to ensure the stability and accuracy of the calculation. This embodiment uses a fixed source term coordinate set to enforce constraints on the iterative intermediate field, that is, it forcibly resets the pixel values ​​located at the fixed source term coordinate set positions in the iterative intermediate field to the corresponding original inversion residual values, preventing the error information of known points from drifting during iteration. After completing the constraints, this embodiment calculates the overall rate of change of the error update matrix and compares it with a preset convergence threshold. If the rate of change is not less than the threshold, this embodiment uses the current iterative intermediate field as the new initial error field and returns to execute the gradient calculation step to enter the next iteration; if the rate of change is less than the threshold, it indicates that the diffusion process has reached a steady state, this embodiment stops iterating and outputs the final iterative intermediate field as a global continuous residual correction surface.

[0037] Furthermore, the specific implementation process of step 600 is as follows: This embodiment performs the final data synthesis operation, that is, by using algebraic superposition to backfill the error information obtained from the physical constraint diffusion solution into the water depth inversion model, thereby generating an initial seabed elevation model with geospatial reference. The construction logic of this model is based on the principle of additive error correction; The specific calculation process is as follows: First, the trend depth value in the relative depth trend surface at the current geospatial plane coordinates is linearly transformed using the linear water depth inversion model determined in the previous steps. This involves multiplying the trend depth value by a transformation gain coefficient and adding a transformation offset coefficient to obtain the predicted water depth surface value at that coordinate. Then, the predicted water depth surface value is arithmetically added to the residual correction value at the corresponding coordinate on the global continuous residual correction surface to finally obtain the initial seafloor elevation model value at that coordinate. In this calculation system, the initial seafloor elevation model is the final output, representing the absolute water depth data field that integrates spectral topographic trends and control point error constraints; the predicted water depth surface is an intermediate variable, representing the theoretical absolute water depth calculated solely through linear regression calibration. The conversion gain coefficient is derived from the least squares fitting method and is used to map the dimensionless relative depth to the meter-level physical depth, for example, a value of 25.5; the trend depth value is derived from the dual-band logarithmic ratio calculation and represents the terrain undulation pattern; the conversion offset coefficient is derived from the regression intercept and is used to correct the reference surface deviation, for example, a value of -1.2; the global continuous residual correction surface is derived from the solution of the non-uniform diffusion equation and refers to the full-field error distribution data after being guided by the texture structure; the geospatial plane coordinates are used to locate the specific pixel position.

[0038] Specifically, the expression for the initial seabed elevation model is: ; in, for Initial seabed elevation model at the location; The predicted water depth surface is obtained by linear regression calibration; This is the conversion gain coefficient; For the relative water depth trend surface in Trend depth value at the location; This is the conversion offset coefficient; This is the global continuous residual correction surface obtained by non-uniform diffusion solution; These are the geographic spatial plane coordinates.

[0039] Furthermore, the specific implementation process of step 700 is as follows: This embodiment first initiates an outlier cleaning procedure based on local statistical features to identify and repair random impulse noise introduced by the data fusion process. This embodiment defines a local sliding window of a preset size on the initial seabed elevation model and controls this window to move pixel-by-pixel across the model matrix, calculating the statistical mean and standard deviation of all valid depth values ​​within the window's coverage area in real time. Based on statistical principles, this embodiment uses the calculated local mean and standard deviation to construct a dynamic outlier threshold, marking pixels whose depth values ​​deviate from the local mean by more than a certain multiple of the standard deviation as isolated noise points. For these marked noise points, this embodiment does not directly remove them as data holes. Instead, it employs inverse distance weighted interpolation, using the valid pixel values ​​in the neighborhood of the noise point as reference nodes. The depth value at this location is numerically reconstructed and replaced according to the inverse relationship of the distance weights, thereby removing spike noise while maintaining the spatial continuity of the seabed topography surface.

[0040] After noise removal, this embodiment further utilizes the Laplacian smoothing operator to perform micromorphological trimming on the denoised elevation model. Since the preceding residual superposition calculation may introduce subtle step-like abrupt changes or unnatural high-frequency oscillations in local areas, this embodiment employs a Laplacian convolution kernel with second-order differential properties to perform iterative convolution processing on the elevation model. Through this iterative smoothing operation, this embodiment appropriately adjusts and converges the depth value of the central pixel based on the depth distribution of neighboring pixels. This effectively weakens and eliminates the step effect on the terrain surface without altering the macroscopic topographic undulation trend, resulting in a more natural and fluid physical characteristic in the microscopic texture of the generated seabed topography.

[0041] Finally, this embodiment performs rigorous spatial boundary constraints and physical attribute compliance checks to ensure the engineering usability of the deliverables. This embodiment uses the binarized water mask of the target image of the water-affected area generated in the previous steps to perform a logical AND operation or spatial cropping on the smoothed elevation model, forcibly removing redundant background data outside the mask area. Simultaneously, this embodiment performs physical logic verification on the depth numerical attribute of each pixel in the model, identifying pixels with positive values ​​as land misclassification pixels or elevation anomalies, and uniformly reassigning these pixels, as well as invalid pixels that were originally non-numerical, to preset invalid data identifier values. After the above series of quality control processes, this embodiment finally outputs a final seabed digital elevation model with clear boundaries, explicit numerical physical meaning, and fully conforming to engineering delivery standards.

[0042] This embodiment also provides a seabed DEM generation system based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, including: The data preprocessing and trend extraction module is used to acquire high-resolution remote sensing images and discrete water depth points of the target sea area, perform radiometric and geometric correction on the high-resolution remote sensing images to obtain target images of the water-bound area, and extract the relative water depth trend surface from the target images of the water-bound area based on the dual-band logarithmic ratio model. The texture structure tensor calculation module is used to calculate the gradient field of the target image of the water-crossing area, obtain the local gradient magnitude and pixel-level structure tensor matrix, and perform feature decomposition on the pixel-level structure tensor matrix to obtain the main feature vector. The diffusion tensor construction module is used to construct an anisotropic diffusion tensor based on the principal eigenvector and the local gradient magnitude. The residual calibration and mesh construction module is used to perform regression calibration on the relative water depth trend surface using the discrete water depth points of the nautical chart within the effective mask range, generate a sparse residual point set, and establish a residual calculation mesh based on the pixel resolution of the target image of the water-bound area; wherein, the effective mask range is a connected region with continuous water spectral characteristics extracted from the target image of the water-bound area using the normalized water index; The non-uniform diffusion solution module is used to construct an initial error field on the residual calculation grid through the sparse residual point set, and use the anisotropic diffusion tensor as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a global continuous residual correction surface. The elevation model synthesis module is used to superimpose the global continuous residual correction surface onto the relative water depth trend surface to obtain an initial seabed elevation model. The quality control and output module is used to remove outliers and smooth edges from the initial seabed elevation model to obtain the final seabed DEM.

[0043] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0044] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, characterized in that, include: High-resolution remote sensing images and discrete water depth points on nautical charts of the target sea area are acquired. Radiometric and geometric corrections are performed on the high-resolution remote sensing images to obtain target images of the water-bound area. Based on a dual-band logarithmic ratio model, a relative water depth trend surface is extracted from the target images of the water-bound area. Gradient field calculation is performed on the target image of the water-contaminated area to obtain the local gradient magnitude and pixel-level structure tensor matrix, and feature decomposition is performed on the pixel-level structure tensor matrix to obtain the main feature vector. An anisotropic diffusion tensor is constructed based on the principal eigenvectors and local gradient magnitudes; Within the effective mask range, the relative water depth trend surface is regressed and calibrated using the discrete water depth points on the nautical chart to generate a sparse residual point set, and a residual calculation grid is established based on the pixel resolution of the target image of the water-intrusion area. The effective mask range is a connected region with continuous water spectral characteristics extracted from the target image of the water-intrusion area using the normalized water index. An initial error field is constructed on the residual calculation grid using the sparse residual point set, and the anisotropic diffusion tensor is used as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a global continuous residual correction surface. The global continuous residual correction surface is superimposed on the relative water depth trend surface to obtain the initial seabed elevation model; Outlier removal and edge smoothing are performed on the initial seabed elevation model to obtain the final seabed DEM.

2. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The process of acquiring high-resolution remote sensing images and discrete depth points from nautical charts of the target sea area, performing radiometric and geometric corrections on the high-resolution remote sensing images to obtain target images of the water-bound area, and extracting relative depth trend surfaces from the target images of the water-bound area based on a dual-band logarithmic ratio model includes: The high-resolution remote sensing image is radiometrically calibrated and atmospherically corrected, converted into a surface reflectance image, and the normalized water index is calculated based on the surface reflectance image. The water-land separation threshold is determined based on the normalized water index, a binary water mask is constructed, and the non-water pixels of the surface reflectance image are removed using the binary water mask to obtain the target image of the water-related area. The blue band reflectance matrix and the green band reflectance matrix are separated from the target image of the water-crossing area; The blue band reflectance matrix and the green band reflectance matrix are subjected to natural logarithmic transformation respectively, and the ratio of the blue and green bands after the logarithmic transformation is calculated to obtain the original relative water depth matrix. The output of the dual-band logarithmic ratio model is the ratio of the natural logarithm of the blue band to the natural logarithm of the green band. The original relative depth matrix is ​​subjected to median filtering for noise reduction to obtain the relative depth trend surface.

3. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The step involves calculating the gradient field of the target image of the water-affected area to obtain the local gradient magnitude and pixel-level structure tensor matrix, and then performing eigenvalue decomposition on the pixel-level structure tensor matrix to obtain the principal feature vector, including: The target image of the water-traveling area is convolved horizontally and vertically using a preset gradient convolution kernel to obtain the horizontal gradient component and the vertical gradient component, respectively. The local gradient magnitude is obtained by performing vector magnitude calculation based on the horizontal gradient component and the vertical gradient component. A gradient autocorrelation matrix is ​​constructed based on the horizontal gradient component and the vertical gradient component, and the gradient autocorrelation matrix is ​​spatially smoothed using a Gaussian convolution kernel to obtain the pixel-level structure tensor matrix. The elements of the gradient autocorrelation matrix include: the square of the horizontal gradient component, the square of the vertical gradient component, and the product of the horizontal gradient component and the vertical gradient component. The pixel-level structural tensor matrix is ​​decomposed into eigenvalues ​​to obtain two eigenvalues ​​and two corresponding eigenvectors. The eigenvector corresponding to the eigenvalue with the larger absolute value is determined as the principal eigenvector.

4. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The construction of the anisotropic diffusion tensor based on the principal eigenvector and the local gradient magnitude includes: Calculate the secondary feature vector based on the primary feature vector; The local gradient magnitude is mapped and calculated using a preset edge stopping function to obtain the normal diffusion coefficient, and a preset tangential diffusion coefficient is set. Construct a normal projection matrix using the principal eigenvectors and the normal diffusion coefficients; A tangential projection matrix is ​​constructed using the secondary eigenvectors and the tangential diffusion coefficients; The anisotropic diffusion tensor is obtained by performing a matrix addition operation between the normal projection matrix and the tangential projection matrix. The expression for the edge stopping function is: ； in, The normal diffusion coefficient is the output of the edge stopping function; The magnitude of the local gradient; The gradient scale parameter is used to characterize the boundary scale of edge intensity; This is a grayscale or single-band reflectance intensity map of the target image in the water-contaminated area. The spatial gradient operator for image intensity; Let be the Euclidean norm of the vector.

5. A method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 4, characterized in that, The step of calculating the secondary feature vector based on the primary feature vector includes: The principal feature vector is analyzed by coordinate analysis to separate the principal horizontal component value and the principal vertical component value; Based on the principle of orthogonality of two-dimensional vectors, the secondary horizontal component value and the secondary vertical component value are determined, wherein the secondary horizontal component value is set as the opposite of the primary vertical component value, and the secondary vertical component value is set as the primary horizontal component value. The secondary horizontal component value and the secondary vertical component value are combined in a vector format to obtain the secondary feature vector.

6. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, Within the effective masking range, the relative depth trend surface is regressed and calibrated using the discrete depth points on the nautical chart, generating a sparse residual point set. A residual calculation grid is then established based on the pixel resolution of the target image of the wading area, including: The discrete water depth points on the nautical chart are traversed, and points whose spatial coordinates are within the effective mask range are selected as effective water depth control points. The actual water depth value corresponding to each effective water depth control point is then extracted. Based on the spatial coordinates of each effective water depth control point, the corresponding trend depth value is extracted from the relative water depth trend surface, and a regression sample pair consisting of multiple sets of trend depth values ​​and actual water depth values ​​is constructed. The least squares method was used to perform linear fitting on the regression sample pairs to determine the conversion gain coefficient and conversion offset coefficient, and a linear water depth inversion model was constructed. The trend depth value is converted into the model predicted depth value using the linear water depth inversion model. The difference between the actual water depth value and the model predicted water depth value is calculated to obtain the inversion residual value. The coordinates of all the effective water depth control points and their corresponding inversion residual values ​​are used to form the sparse residual point set. Based on the sparse residual point set, the number of rows, columns and geospatial resolution of the target image of the water-crossing area are obtained, a zero-value matrix of the same dimension as the target image of the water-crossing area is constructed, and the zero-value matrix is ​​marked as the residual calculation grid. The expression for the linear water depth inversion model is: ； in, coordinates The model predicts the water depth at that location; For the relative water depth trend surface in Trend depth value at the location; The conversion gain coefficient is used to map the trend depth scale to the water depth scale; The conversion offset coefficient is used for reference translation; These are the geographic spatial plane coordinates.

7. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The process involves constructing an initial error field on the residual calculation grid using the sparse residual point set, and using the anisotropic diffusion tensor as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a globally continuous residual correction surface, including: Each inversion residual value in the sparse residual point set is mapped to the corresponding pixel position of the residual calculation grid, and the unmapped pixel positions are initialized to zero to obtain the initial error field. At the same time, the pixel positions mapped with the inversion residual values ​​are marked as a fixed source term coordinate set. Gradient calculation is performed on the initial error field to obtain the error gradient field, and matrix multiplication is performed between the anisotropic diffusion tensor and the error gradient field to obtain the diffusion flux vector field. Divergence calculation is performed on the diffusion flux vector field to obtain the error update matrix, and the initial error field is numerically accumulated and updated according to the preset time step and the error update matrix to obtain the iterative intermediate field; The iterative intermediate field is subjected to forced constraints using the fixed source term coordinate set; Determine whether the overall rate of change of the error update matrix is ​​less than a preset convergence threshold; if it is not less than a preset convergence threshold, then use the iterative intermediate field as the new initial error field and return to the step of performing the gradient operation; if it is less than a preset threshold, then output the iterative intermediate field as the global continuous residual correction surface.

8. The method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The expression for the initial seabed elevation model is: ； in, for Initial seabed elevation model at the location; The predicted water depth surface is obtained by linear regression calibration; This is the conversion gain coefficient; For the relative water depth trend surface in Trend depth value at the location; This is the conversion offset coefficient; This is the global continuous residual correction surface obtained by non-uniform diffusion solution; These are the geographic spatial plane coordinates.

9. A method for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces according to claim 1, characterized in that, The process of removing outliers and smoothing edges from the initial seabed elevation model to obtain the final seabed DEM includes: The initial seabed elevation model is subjected to statistical analysis based on local neighborhood to calculate the mean depth and standard deviation of each pixel within a preset neighborhood window; An anomaly detection threshold is set based on the mean depth and the standard deviation. Pixels in the initial seabed elevation model that exceed the anomaly detection threshold are marked as noise points. The noise points are then numerically reconstructed using inverse distance weighted interpolation to obtain the denoised elevation model. The denoised elevation model is iteratively smoothed using the Laplace smoothing operator to eliminate the step effect on the terrain surface, resulting in a smoothed elevation model. The smoothed elevation model is cropped using a binarized water mask of the target image of the water-involved area. Misidentified land pixels with positive depth values ​​and invalid pixels with non-numerical depth values ​​are uniformly assigned invalid data identifiers to obtain the final seabed DEM that meets the engineering delivery standards.

10. A system for generating a seabed DEM based on large-scale nautical chart depth points and high-resolution remote sensing surfaces, characterized in that, include: The data preprocessing and trend extraction module is used to acquire high-resolution remote sensing images and discrete water depth points of the target sea area, perform radiometric and geometric correction on the high-resolution remote sensing images to obtain target images of the water-bound area, and extract the relative water depth trend surface from the target images of the water-bound area based on the dual-band logarithmic ratio model. The texture structure tensor calculation module is used to calculate the gradient field of the target image of the water-crossing area, obtain the local gradient magnitude and pixel-level structure tensor matrix, and perform feature decomposition on the pixel-level structure tensor matrix to obtain the main feature vector. The diffusion tensor construction module is used to construct an anisotropic diffusion tensor based on the principal eigenvector and the local gradient magnitude. The residual calibration and mesh construction module is used to perform regression calibration on the relative water depth trend surface using the discrete water depth points of the nautical chart within the effective mask range, generate a sparse residual point set, and establish a residual calculation mesh based on the pixel resolution of the target image of the water-bound area; wherein, the effective mask range is a connected region with continuous water spectral characteristics extracted from the target image of the water-bound area using the normalized water index; The non-uniform diffusion solution module is used to construct an initial error field on the residual calculation grid through the sparse residual point set, and use the anisotropic diffusion tensor as the diffusion coefficient tensor to perform non-uniform diffusion calculation on the initial error field to obtain a global continuous residual correction surface. The elevation model synthesis module is used to superimpose the global continuous residual correction surface onto the relative water depth trend surface to obtain an initial seabed elevation model. The quality control and output module is used to remove outliers and smooth edges from the initial seabed elevation model to obtain the final seabed DEM.

Images