Large scene water depth model reconstruction method based on sparse single-beam sounding data constraint

Abstract

The application discloses a kind of based on sparse single-beam sounding data constraint large scene water depth model reconstruction method, including the following steps: single-beam sounding point cloud is rasterized and cut as low-quality data, with complete DEM as high-quality data, with topographic feature line as model priori knowledge, the data set required for model training and test is made;The process of bidirectional diffusion is modeled in the form of Gaussian distribution;Condition denoising network is constructed, and prior information is encoded and injected into denoising main network by condition information encoder;The trained model is used for cyclic sampling, the sampling result is carried out single-step forward diffusion, and the single-step sampling step is repeated multiple times to eliminate the semantic conflict between measured data and sampling result;The image block after reconstruction is spliced.The present application can reconstruct complete high-precision seabed topographic model according to sparse discrete single-beam sounding data, and can play a practical value when global sounding data coverage is insufficient.

Description

TECHNICAL FIELD

[0001] The present application relates to seabed bathymetric data processing, digital elevation model (DEM) generation and terrain reconstruction technology based on deep learning, and in particular to a large-scene water depth model reconstruction method based on sparse single-beam bathymetric data constraint. BACKGROUND

[0002] Seafloor digital elevation model is a scientific expression of the shape of the seafloor surface, which is crucial in many fields such as global plate tectonics, marine mineral survey, benthic biological survey, and marine geological activity. However, electromagnetic wave measurement medium has the disadvantage of strong attenuation in seawater, making it difficult for traditional surveying and mapping techniques to work in deep sea topography measurement. Considering the low efficiency of ship-borne sonar sounding and long implementation period, bathymetric blank areas exist widely in the world, and this problem is difficult to improve significantly in a short period of time. In order to accelerate the construction of global high-resolution DEM, GEBCO launched the "Seabed 2030" project, whose goal is to compile a seamless global seafloor digital map by 2030. Single-beam sounding technology originated early in the history of human exploration of the ocean, and has accumulated a large amount of data, which is an important source of measured data for global seafloor topographic modeling. However, these data show random and sparse spatial distribution characteristics, which limit their continuous expression of seafloor topography and geomorphology.

[0003] Due to the original measurement data, the problem of DEM producing voids is inevitable. Overcoming the sparsity of data and accurately restoring the true terrain from incomplete measurement data has a long history. Initially, researchers used spatial interpolation techniques to estimate the elevation of unknown areas, such as inverse distance weighting, Kriging, and triangulated irregular network methods. The premise of using interpolation methods to estimate the missing parts in DEM is to assume that the terrain changes are continuous, and the elevation of unknown points is highly correlated with the neighborhood points. However, traditional interpolation methods have limited ability to capture feature changes and are only effective in small voids or flat terrain. In complex terrain or areas with large data gaps, it often causes the terrain to be too smooth and locally distorted.

[0004] In recent years, deep learning technology has made breakthroughs and gradually developed into a new paradigm for solving image inpainting problems, providing a promising direction for DEM inpainting. Especially the generative adversarial network (GAN) can generate clear and complete semantic structures from sparse data, improving the limitations of interpolation methods in handling complex terrain DEM inpainting. In particular, the self-attention mechanism is introduced into GAN, which enables the network to better learn the relevance of the context and achieve optimal results in various terrain backgrounds. However, this unconditional repair mode has strong randomness, and researchers try to incorporate more prior knowledge into the model, mainly including two categories: one is quantitative information, including remote sensing images, terrain shadows, slopes, and relief textures; the second is qualitative information, mainly including ridge lines, valley lines, and other terrain feature lines. These prior terrain knowledge effectively improves the consistency of the DEM inpainting results. Although GAN can learn non-linear mapping relationships in complex visual tasks, the mode collapse problem is difficult to completely solve, resulting in difficulty in training GAN models, low scalability, and the generated results may contain unreasonable artifacts. Therefore, the work of DEM reconstruction from sparse bathymetric data sources still needs to be further improved. SUMMARY

[0005] In order to overcome the shortcomings of the prior art, the purpose of the present application is to provide a large scene water depth model reconstruction method based on sparse single-beam bathymetric data constraint. For the first time, sparse single-beam bathymetric data is used as prior constraint, terrain feature lines are used as conditional input, and the cyclic sampling elimination spatial inconsistency mechanism is combined to achieve high-precision, seamless large scene reconstruction. Through testing, the present application can reconstruct high-precision DEM from only 5% of the water depth points, and plays a role in the process of seafloor topography measurement and modeling.

[0006] In order to achieve the above purpose, the present application adopts the following technical solutions:

[0007] A large scene water depth model reconstruction method based on sparse single-beam bathymetric data constraint, comprising the following steps:

[0008] (1) Perform rasterization processing on the single-beam bathymetric point cloud to obtain low-quality data (low quality, denoted as LQ) containing bathymetric blank areas; take the complete digital elevation model DEM as high-quality data (high quality, denoted as HQ); extract the terrain feature lines to construct the ridge line, valley line and transition zone terrain feature line layer C; spatially register LQ, HQ and C, and divide them into image blocks with a size of 120x120 pixels and an adjacent overlap rate of 50%;

[0009] (2) Use a conditional diffusion probability model, and the forward diffusion process is defined as: at time step t, add Gaussian noise to HQ to obtain (1),

[0010] where ε~N(0, 1), β t is a preset noise variance sequence, t = 1, 2, …, T; the backward diffusion process is predicted by the neural network ε θ , and the reconstructed x t is obtained from x t−1 ;

[0011] (3) A conditional denoising network is constructed, whose backbone is a four-layer U-shaped network structure (denoted as UNet), each layer containing two dense residual modules and a scale change layer; a Transformer module is embedded at the bottom of the UNet to model long-range dependencies; independent encoders for LQ and C are constructed respectively, each being a three-layer convolution-residual-downsampling structure; the outputs of the two encoders are concatenated in the channel dimension and then injected into the bottom Transformer of the UNet to guide noise prediction;

[0012] (4) In the inference stage, forward diffusion is performed on LQ to obtain x t ; at each time step t, the known areas in x t are kept as the original LQ values based on the mask (a binary image recording the positions of the measured points), and the unknown areas are predicted by the conditional denoising network to obtain ε θ , which is then updated by backward diffusion to generate x t−1 ; this process is repeated N times or until the relative mean square error between the known areas and the predicted areas is less than a preset threshold δ to eliminate spatial inconsistencies;

[0013] (5) All reconstructed image blocks are horizontally concatenated and then vertically concatenated; in the overlapping area (width L pixels) of adjacent image blocks, a weighted blending function is used to calculate the elevation value: (12), where k is the pixel distance of the point from the left edge, D A , D B are the elevation values at the corresponding positions of the adjacent blocks, and the complete DEM of the large scene is finally generated.

[0014] The method, in step (1), the resampling resolution r is 5m, 10m or 30m; the topographic feature line layer C is generated by extracting the ridge line and the valley line from the General Bathymetric Chart of the Oceans (GEBCO) through the hydrological analysis tool of ArcGIS.

[0015] The method, in step (1), the topographic feature line layer C is constructed by ridge line pixel value 1, valley line pixel value 0 and transition zone pixel value 0.5.

[0016] The method, in step (2), β t The cosine scheduling strategy is adopted: (2),

[0017] Wherein T =1000, β min=0.0001, β max=0.02.

[0018] The method, in step (3), the Transform module is an 8-head self-attention structure, and the hidden layer dimension is 256.

[0019] The method, in step (4), the mask is generated by the non-empty value area in the LQ, that is, if the depth value of a point exists, the mask=1; otherwise, the mask=0.

[0020] The method, in step (4), the relative mean square error threshold δ is 0.05 meters, or the maximum iteration number N is 30 times.

[0021] An apparatus for implementing the method comprises a processor, a memory and an input / output interface; the processor is configured to execute all the steps; the memory is used to store the LQ, the HQ, the C, the model parameters and the intermediate results; and the input / output interface is used to receive the single-beam sounding point cloud and output the reconstructed large scene DEM.

[0022] A computer readable storage medium having a computer program stored thereon, when the program is executed by a processor, the method is implemented.

[0023] The beneficial technical effects brought by the present application are:

[0024] This invention discloses a method for reconstructing large-scene seabed depth models based on sparse single-beam bathymetry data constraints. Compared with existing technologies, this invention predicts seabed depth values ​​in blank areas of single-beam bathymetry based on probabilistic generation theory. Simultaneously, a conditional encoder is designed to enhance the role of prior knowledge in the DEM reconstruction process, effectively limiting the diversity of generated results and achieving more accurate DEM reconstruction results. Furthermore, this invention achieves seabed DEM reconstruction in large-scene scenarios by processing edge errors. This invention will further promote the construction of global seabed DEMs based on measured data, playing an important role in marine engineering construction, protection of maritime rights, seabed resource development, and marine surveying. Attached Figure Description

[0025] Figure 1 This is the flowchart of the large-scene water depth model reconstruction method based on sparse single-beam bathymetry data constraints in this invention.

[0026] Figure 2 This describes the bidirectional diffusion modeling process in this invention.

[0027] Figure 3 This is a diagram of the conditional noise reduction network structure in this invention.

[0028] Figure 4 This is the cyclic sampling process in this invention.

[0029] Figure 5 This is the principle of edge joint error handling in this invention. Detailed Implementation

[0030] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0031] This study uses both simulated and measured single-beam bathymetry (LQ) data. The simulated data involves randomizing the complete DEM using simulated random tracklines, while the measured data uses publicly available global single-beam bathymetry point clouds. The 2024 version of the general ocean depth map GEBCO_2024 is used as the HQ, and topographic feature lines are generated using ArcGIS hydrological analysis tools. The overall workflow is as follows: Figure 1The method comprises the following steps: step 1, rasterizing and cutting single-beam sounding point cloud as low-quality data LQ, taking a complete DEM as high-quality data HQ, taking a topographic feature line as model prior knowledge, and preparing a data set required for model training and testing; step 2, modeling a bidirectional diffusion process in the form of a Gaussian distribution, wherein the forward diffusion is gradually adding Gaussian noise in the complete DEM to expand the boundary of the distribution, and the reverse diffusion is constantly shrinking the boundary of the distribution by predicted noise to finally complete the prediction of the target domain; step 3, constructing a conditional denoising network, encoding prior information through a conditional information encoder and injecting the denoising backbone network to realize the guidance of the diffusion process and complete the prediction of the noise added in the forward diffusion process; step 4, using the trained model to perform cycle sampling, implementing forward diffusion on LQ at a specific time step, implementing reverse diffusion on the noise image, and fusing the two as a single sampling result, performing single-step forward diffusion on the sampling result, and repeating the single-step sampling step multiple times to eliminate the semantic conflict between the measured data and the sampling result; and step 5, constructing a blending function to process the edge error of the repaired small-size DEM, and splicing the small-size repaired result into a complete DEM in a large scene.

[0032] Step 1: rasterizing and cutting single-beam sounding point cloud as low-quality data LQ, taking a complete DEM as high-quality data HQ, and taking a topographic feature line as model prior knowledge to prepare a data set required for model training and testing, and the method is as follows:

[0033] Step 1.1: first, rasterize the discrete single-beam sounding point cloud, at this time, the data contains a large number of sounding blank areas, denoted as low-quality data LQ, at this time, a resampling resolution needs to be set according to the point cloud data distribution r , and the resolution will also determine the resolution of the final reconstructed DEM.

[0034] Step 1.2: according to the data reflecting the topographic relief changes such as global public bathymetric model or gravity satellite bathymetric inversion model, the ridge line and valley line are made, and the ridge line and valley line are merged into the topographic feature line layer, in the layer, the pixel value of the ridge line position is 1, the pixel value of the valley line position is 0, and the pixel value of the remaining position is 0.5.

[0035] Step 1.3: taking a complete reference DEM as high-quality data HQ, registering LQ, HQ and topographic feature line according to spatial coordinates, and cutting them into 120*120 small images, and reserving 50% overlap rate between adjacent small images to facilitate subsequent processing of edge error.

[0036] Step 2: as Figure 2The process of bidirectional diffusion is modeled in the form of Gaussian distribution, where the forward diffusion is to gradually add Gaussian noise in the complete DEM, expanding the boundary of its distribution, and the reverse diffusion is to constantly shrink the boundary of the distribution by predicted noise, finally completing the prediction of the target domain, the method is as follows:

[0037] Step 2.1: The whole diffusion process is designed as a Markov process, that is, t The state at time t x t Only related to t The state at time t-1 x t-1 , by controlling x t-1 Adding random Gaussian noise in x t , we can get q ( x t | x t-1 ) obeys Gaussian distribution, define the single-step forward diffusion process as:

[0038] (3),

[0039] In the formula, β t is the hyperparameter to control the variance of added noise, and T is the total length of Markov chain, that is, the total number of diffusion steps. According to formula (3), we can further deduce the distribution formula of x 0 to the state x t at any time, that is, the multi-step forward diffusion process, that is:

[0040] (4),

[0041] In the formula, .

[0042] Step 2.2: The reverse diffusion process is to predict the state x t at time t-1 according to the state x t-1 at time t, and this process is recorded as P θ ( x t | x t-1 ), P θ ( x t | x t-1 ) is also Gaussian distribution, that is:

[0043] (5),

[0044] In the formula, ε θ denotes the noise predicted by the neural network, and denotes P θ ( x t | x t-1 ) is the variance, and its calculation formula is:

[0045] (6),

[0046] Step 3: As shown in Figure 3 , a conditional denoising network is constructed, the prior information is encoded by a conditional information encoder and injected into the denoising backbone network, the guidance of the diffusion process is realized, and the prediction of the added noise in the forward diffusion process is completed, the method is as follows:

[0047] Step 3.1: The backbone of the denoising network is designed as a 4-layer UNet structure, as shown in part (a) of Figure 3 , each layer contains two dense residual modules and a scale change layer, and the left side of the UNet is a down-sampling layer and the right side is an up-sampling layer. Since the time complexity of self-attention learning is positively correlated with the square of the image size, in order to improve the efficiency of model training and inference, a Transformer module is set at the bottom of the UNet space for self-attention learning on small size images.

[0048] Step 3.2: The conditional encoder is used to encode the LQ and the terrain feature line. The conditional encoder is designed as a three-layer structure, each layer contains a residual module and a down-sampling module. The LQ and the terrain feature line encoders are constructed respectively, as shown in part (b) of Figure 3 , they are independently encoded, then the encoded information is spliced in the channel direction, and then injected into the Transformer module at the bottom of the backbone network, and the noise prediction learning is performed together with the backbone network.

[0049] Step 3.3: The model is trained according to the training set constructed in step 1, and the objective function is:

[0050] (7),

[0051] In the formula, θ denotes the learnable parameters in the neural network, x t denotes the noise added in the forward diffusion process, ε θ x t ,c,t,r represents the noise prediction result, c represents the ridge line and the valley line, t represents the current time step, r represents the preset resolution during resampling.

[0052] Step 4: using the trained model for cycle sampling, as shown in Figure 4 , the LQ is forward diffused at a specific time step, the noisy image is backward diffused, and the two are fused as a single sampling result. The sampling result is single-step forward diffused, and the single-step sampling step is repeated multiple times to eliminate the semantic conflict between the measured data and the sampling result, and the method is as follows:

[0053] Step 4.1: forward diffuse the LQ according to formula (4) to calculate x t-1 The part of the measured data is denoted as , which can be directly calculated from the known data, and the calculation formula is:

[0054] (8),

[0055] In the formula, ε ~N(0,1), is the result obtained by sampling from the standard Gaussian distribution, and mask is an image recording the position of the measured water depth. The pixel value of the position with measured water depth is 1, and the pixel value of the position without measured water depth is 0.

[0056] Step 4.2: calculate the backward diffusion of x t according to formula (5) to predict x t-1 The part of the measured data in is denoted as

[0057] (9),

[0058] In the formula, ε'~N(0,I) is the result obtained by sampling from the standard Gaussian distribution.

[0059] Step 4.3: integrate the results obtained in steps 4.1 and 4.2 into an image to obtain a single sampling result x t-1 , that is:

[0060] (10),

[0061] Since​x t-1 The two parts are obtained using different methods, and there may be semantic conflicts between them. In this case, the content obtained by formula (10) x t-1 Diffusion is carried out one step forward to obtain x t The calculation formula is as follows:

[0062] (11),

[0063] Step 4.4: In the new round of sampling, the result obtained from formula (11) x t Substitute the values ​​back into formula (5) for reverse diffusion, and repeat steps 4.2-4.3 until... x t-1 The semantic inconsistency disappears. For example, when the relative mean square error between the known region (mask=1) and the predicted region is less than 0.05 (i.e. δ=0.05m) in two consecutive iterations, or when the number of iterations reaches 30, the cyclic sampling is terminated.

[0064] Step 5: Construct a blending function to handle the edge joint errors of the repaired small-size DEM. The overall process is as follows: Figure 5 As shown, the method for stitching small-sized repair results into a complete DEM for a large scene is as follows:

[0065] Step 6.1: In step 1.3, a 50% overlap rate was reserved when cropping the small-sized DEM. The elevation of the overlapping area was predicted twice. A blending function was designed based on trigonometric functions for edge joining, and the weighted average of the elevations within the overlapping area was calculated. Let the two adjacent small-sized restoration results be images... A and images B The formula for calculating the elevation of a point within the mixed region is as follows:

[0066] (12),

[0067] In the formula, D This indicates the final elevation value of the point. k This indicates the distance between the point and the left edge of the mixing area. L Indicates the width of the blending region. D A Representing an image A The elevation value of that point in the data. D B Representing an image BThe point elevation value in the point. The mixed function is normalized weighting, and the weight sum is constant 1 under any k∈[0, L]. The function realizes smooth transition, avoids the 'platform effect' of linear interpolation in the center of the overlapping area, and improves the continuity of steep slope terrain.

[0068] Step 6.2: transverse splicing is performed on all small size repair results according to step 6.1.

[0069] Step 6.3: all transverse splicing results are vertically spliced according to step 6.1, and finally the complete DEM under the large scene is obtained. The method of the application is used for precision evaluation in multiple areas, and the results are summarized in Table 1. Compared with the inverse distance weighted interpolation, Kriging interpolation and irregular triangle network interpolation method, the application achieves the optimal effect in all areas, indicating that the method of the application is highly close to the real terrain.

[0070] Table 1

[0071]

[0072] The above-described embodiments only express several embodiments of the present application, and the description is more specific and detailed, but it cannot be understood as a limitation on the scope of the application. It should be noted that for ordinary skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which all belong to the protection scope of the present application. The protection scope of the present application is given by the appended claims and any equivalent technical solutions thereof.

Claims

1. A method for reconstructing a large-scene water depth model based on sparse single-beam bathymetry data constraints, characterized in that, Includes the following steps: (1) Rasterize the single-beam bathymetry point cloud to obtain low-quality data LQ containing the bathymetry blank area; use the complete digital elevation model DEM as high-quality data HQ; extract the terrain feature lines and construct the terrain feature line layer C of ridge line, valley line and transition area; spatially register LQ, HQ and C and divide them into image blocks with a size of 120×120 pixels and an adjacent overlap rate of 50%; (2) Using the conditional diffusion probability model, the forward diffusion process is defined as follows: at time step t, Gaussian noise is gradually added to HQ to obtain , Where ε ~ N(0,1), β t A pre-defined noise variance sequence is given, t=1,2,…,T; the back diffusion process is predicted by a neural network. ε θ To achieve from x t reconstruction x t−1 ; (3) Construct a conditional noise reduction network with a four-layer U-shaped network structure UNet as its backbone. Each layer contains two dense residual modules and one scale-changing layer. Embed a Transformer module at the bottom of UNet to model long-range dependencies. Construct independent encoders for LQ and C respectively, each with a three-layer convolution-residual-downsampling structure. Concatenate the outputs of the two encoders in the channel dimension and inject them into the bottom Transformer of UNet to guide noise prediction. (4) During the inference phase, forward diffusion is performed on LQ to obtain x t At each time step t, a binary image recording the measured point locations based on the mask will be generated. x t The known regions retain their original LQ values, while the unknown regions are predicted using a conditional denoising network. ε θ And perform reverse diffusion update to generate x t−1 ; Repeat this process N times or until the relative mean square error between the known region and the predicted region is less than a preset threshold δ, in order to eliminate spatial inconsistency. These include: Step 4.1: Perform forward diffusion on LQ and calculate... x t-1 The portion of the measured data is denoted as The calculation formula is: (8), where, The mask is an image that records the location of the measured water depth. The pixel value is 1 for the location of the measured water depth and 0 for the location of the location of the measured water depth. Step 4.2: For x t Perform reverse diffusion calculations to predict x t-1 The blank area in the medium depth sounding is denoted as The calculation formula is: （9）； In the formula (6); Step 4.3: Combine the results obtained in Steps 4.1 and 4.2 into a single image to obtain the single sampling result. x t-1 ,Right now: （10）； (5) All reconstructed image blocks are stitched horizontally and then vertically; in the overlapping area of ​​adjacent image blocks, with a width of L pixels, the elevation value is calculated using a weighted mixing function: , in k This is the pixel distance from the current point to the left edge. D A , D B The elevation values ​​of the corresponding locations of adjacent blocks are used to generate a complete large-scene DEM.

2. The method as described in claim 1, characterized in that: In step (1), the resampling resolution r is 5m, 10m or 30m; the topographic feature line layer C is generated by extracting ridge lines and valley lines from the general ocean depth map GEBCO using ArcGIS hydrological analysis tools.

3. The method as described in claim 1, characterized in that: In step (1), a terrain feature line layer C is constructed with a ridgeline pixel value of 1, a valley line pixel value of 0, and a transition area pixel value of 0.

5.

4. The method as described in claim 1, characterized in that: In step (2), β t Cosine scheduling strategy adopted: , in T =1000, β min =0.0001, β max =0.

02.

5. The method as described in claim 1, characterized in that: In step (3), the Transformer module is an 8-head self-attention structure with a hidden layer dimension of 256.

6. An apparatus for implementing the method according to any one of claims 1-5, characterized in that: The system includes a processor, a memory, and an input / output interface; the processor is configured to perform all the steps of the method described in any one of claims 1-5; the memory is used to store LQ, HQ, C, model parameters, and intermediate results; the input / output interface is used to receive single-beam bathymetry point clouds and output the reconstructed large-scene DEM.

7. A computer-readable storage medium having a computer program stored thereon that, when executed by a processor, implements the method of any one of claims 1-5.

Images