Intelligent recognition and deformation prediction method for submarine landslide based on time-series multi-beam data

Abstract

The application provides a submarine landslide intelligent identification and deformation prediction method based on time sequence multi-beam data. The method takes 2.5 m grid DEM obtained by a ship-borne multi-beam sounding system as input, and automatically extracts single-stage submarine landslide rear edge vector through GST coherence weighting, 3D weighted principal component analysis, three-dimensional Kuwahara filtering and region growing algorithm in turn; then, confidence, trend and geometric type are obtained simultaneously by using feature enhancement vectorization, three-dimensional mapping and attribute calculation, so that fine identification of submarine landslide rear edge is realized. Spatial registration based on features is carried out on different period rear edge lines, displacement, trend and curvature difference are calculated point by point, and finally through time and space sliding window weighted extrapolation, historical difference features are superimposed to the latest rear edge line, the next period prediction rear edge vector is automatically generated, and the automatic identification and deformation prediction of multi-period submarine landslide are completed. The method has the advantages of full automation, few parameters, strong noise resistance, short cycle, integrated identification and prediction and the like.

Description

Technical Field

[0001] This invention relates to the field of submarine landslide identification and deformation prediction technology, and more specifically, to a method for intelligent identification and deformation prediction of submarine landslides based on temporal multibeam data. This invention automatically extracts the trailing edge vector of submarine landslides using multi-phase multibeam GST, 3D-WPCA, Kuwahara, and region growing algorithms, and achieves dynamic monitoring and intelligent prediction through feature registration and spatiotemporal sliding window extrapolation. Background Technology

[0002] Submarine landslides are gravity flow processes that involve the transport of loose seabed sediments along continental slopes in whole or in discrete blocks, triggered by factors such as earthquakes, tsunamis, rapid sedimentation, volcanic activity, hydrate decomposition, or tidal erosion. These processes include sliding, collapse, and debris flows. This hazard is widespread along both passive and active continental margins, with individual landslides ranging in size from a few square meters to thousands of square kilometers, posing a significant threat to drilling platforms, submarine pipelines, fiber optic cables, and port infrastructure.

[0003] Traditional landslide surveys rely on single multibeam bathymetry and manual interpretation of 2D or 3D seismic profiles, which has three major bottlenecks: high data acquisition costs, as it is difficult to repeatedly collect high-precision multibeam and seismic data in the same area; long interpretation cycles, requiring senior geophysicists to compare each survey line, resulting in low efficiency; and significant information loss after the identification results are normalized, coupled with a simplistic principle, low data utilization, and difficulty in continuously tracking the dynamic changes of landslide bodies, leading to delayed risk warnings.

[0004] With the advancement of the national marine strategy, near-shore engineering is showing a trend towards grid-based and intensive development, posing new demands for real-time, comprehensive, and dynamic monitoring of submarine landslides. In recent years, shipborne multibeam bathymetry systems have enabled repeated measurements of the same sea area annually or even quarterly, forming time-series DEM data, providing a data foundation for multi-period comparisons. However, the automatic processing of massive amounts of water depth data, accurate conformal registration between different periods, and reliable extraction of subtle topographic changes remain pressing technical challenges that need to be addressed.

[0005] While deep learning has achieved remarkable results in image recognition, the scarcity of submarine landslide samples and the high cost of annotation make it difficult to implement purely data-driven models. To address this, this invention proposes an intelligent identification and deformation prediction method for submarine landslides based on temporal multibeam data. Using GST coherence-weighted three-dimensional weighted principal component analysis (3D-WPCA) as its core, it automatically extracts the trailing edge vector of the current landslide without requiring any slope threshold. Then, improved ICP is used to complete multi-period vector spatial registration, and the predicted trailing edge line for the next period is automatically generated through temporal sliding window TSW weighted extrapolation. The entire process eliminates the need for manual delineation of the slip zone, enabling intelligent identification and deformation prediction of multi-period, multibeam submarine landslides, providing efficient and low-cost dynamic data support for the full lifecycle safety management of nearshore engineering projects. Summary of the Invention

[0006] To overcome the shortcomings of traditional techniques such as low computational accuracy and insufficient data utilization, this invention provides a time-series multibeam bathyspherical landslide identification and deformation prediction method that integrates single-phase intelligent extraction, multi-phase spatial registration, and spatiotemporal model prediction. This method uses a 2.5 m grid DEM sequence repeatedly acquired by conventional shipborne multibeam bathyspherical systems as the data source. In a single phase, it sequentially performs GST coherence weighting, 3D-WPCA, Kuwahara filtering, region growing, vectorization extraction, and attribute calculation. In multiple phases, spatial registration is completed through improved ICP, and historical difference features are weighted and extrapolated using a spatiotemporal sliding window (TSW) to automatically generate the predicted trailing edge line for the next phase. This achieves intelligent identification and prediction of the entire process of submarine landslides without manual delineation.

[0007] This invention is achieved through the following technical solution: a method for intelligent identification and deformation prediction of submarine landslides based on time-series multibeam data, characterized by the following steps: Step S1: Acquisition and preprocessing of raw data in multiple periods T0, T1, and T2. The target sea area is measured using the same shipborne multibeam system to obtain raw data. After acoustic ray tracking, tide level correction, and noise reduction, a 2.5 m × 2.5 m regular grid DEM is output, namely DEM_T0, DEM_T1, and DEM_T2. Step S2: Calculate the coherence of the DEM for each single period using the gradient structure tensor GST to obtain the reliability weight of each grid point. Step S3: Using the coherence from Step S2 as weights, perform 3D weighted principal component analysis (3D-WPCA) within a 5×5×5 window to obtain the gradient and dip angle, generating an initial landslide trailing edge probability volume; apply 3D Kuwahara filtering to the probability volume, preserving edges and suppressing noise according to the principle of minimizing the variance of (2n+1)×(2n+1) sub-windows, to obtain an enhanced landslide trailing edge probability map; automatically select seed points on the probability map, execute a region growing algorithm, and aggregate similar probability pixels into continuous linear objects to complete the initial extraction of the landslide trailing edge; Step S4: Reduce the dimensionality of the 3D probabilistic volume to a 2D plane, and use a feature enhancement algorithm to normalize and enhance the probabilistic map, making the pixel values ​​at the landslide trailing edge polarized; after skeletonization, broken line connection and spur removal, the enhanced binary image is processed to obtain a single-pixel wide continuous landslide trailing edge line; vectorization tracking is performed along the skeleton line to generate independent vector line segments, and the segments are reprojected onto the original 3D DEM space according to the spatial coordinates to realize the vectorization extraction and 3D mapping of each submarine landslide trailing edge; Step S5: Automatically calculate the confidence level, direction and geometric type attributes of each vectorized landslide trailing edge vector to achieve refined extraction of landslide trailing edge; Step S6: Using the trailing edge vector of the first landslide as the main reference, the improved ICP algorithm is used to spatially register the trailing edge lines of subsequent phases. The registration cost function considers both Euclidean distance minimization and strike consistency constraints. After registration, the displacement vector between adjacent phases is calculated point by point along the normal direction of the trailing edge line at a grid resolution of 2.5 m, and the strike change and curvature change are extracted simultaneously. Step S7: Introduce a spatiotemporal sliding window (TSW) to perform a weighted average of the difference features for each single period. The weights are determined by the confidence level and the time distance from the prediction period to obtain the predicted difference features. The predicted difference features are superimposed on the latest period trailing edge line, and the trailing edge line of the prediction period is generated by extrapolating point by point to realize the prediction of submarine landslide deformation from time-series multibeam data.

[0008] As a preferred option, the GST coherence calculation in step S2 specifically includes the following steps: Step S21: Calculate ∂h / ∂x, ∂h / ∂y, and ∂h / ∂z pixel by pixel for DEM_Ti to construct the gradient vector ∇h; Step S22: Accumulate the outer product in the 3×3×3 neighborhood to obtain the gradient structure tensor: GST = Σ (∇h)(∇h)^T; Step S23: Find the eigenvalues ​​λ1≥λ2≥λ3, and calculate coherence_GST = (λ1−λ3) / λ1, taking values ​​from 0 to 1. The closer to 1, the more uneven the local surface is, which is used as the reliability weight for subsequent 3D-WPCA.

[0009] As a preferred option, the three-dimensional weighted principal component analysis in step S3 specifically includes the following steps: Step S311: Run WPCA point-by-point on DEM_Ti using a 5×5×5 voxel sliding window (actual area ≈ 12.5 m × 12.5 m × 5 m), where C = Σ w · (∇h)(∇h)^T, weight w = exp(−d / d0) · coherence_GST, d0=2 voxels; Step S312: Extract the minimum eigenvalue λ_min of C, which is defined as the roughness index U; U greater than 0.65 is considered as a potential trailing edge voxel; Step S313: The output is a three-dimensional probability volume Probability_Ti, with a value of 0-1. The larger the value, the more likely it is to be the trailing edge of the landslide.

[0010] Furthermore, the three-dimensional Kuwahara filtering in step S3 preserves edges and suppresses noise, specifically including the following steps: Step S321: Divide Probability_Ti into 8 overlapping 3×3×3 sub-blocks; Step S322: Calculate the mean m_k and standard deviation σ_k of each sub-block, and retain the mean of the sub-block corresponding to σ_min as the filter output; Step S323: Preserve the true trailing edge abrupt change signal and suppress side noise.

[0011] Furthermore, the region growing clustering in step S3 specifically includes the following steps: Step S331: Automatic seed selection: Find connected cores with local maxima and U>0.75 within the filtered U-body; Step S332: Use the similarity criterion |U_i – U_seed|<0.1 to make a judgment; only those with a continuous number of voxels ≥200 m² are retained. Step S333: The output is a raster cluster Cluster_Ti, with each cluster assigned a unique value.

[0012] As a preferred option, step S4 specifically includes the following steps: Step S41: Project Cluster_Ti onto a two-dimensional plane with the highest probability along the depth direction; Step S42: Use grayscale remapping I′ = (I–I_min)^γ / (I_max–I_min)^γ, γ=0.5 to achieve pixel polarization; Step S43: First, binarize the image to a threshold of 0.5 to obtain a binary image B. After refining image B into a single-pixel wide skeleton, automatically connect the broken ends with a Hausdorff distance of less than 3 pixels and a direction difference of less than 30°. Automatically remove burrs with a length of less than 5 pixels. Only retain single trailing edge lines with a length of 40 m or more.

[0013] As a preferred option, step S5 specifically includes the following steps: Step S51: Map the vector of the trailing edge of a single landslide back to the original 3D DEM. Extract the GST coherence, 3D-WPCA gradient magnitude, and Kuwahara variance within the buffer zone. After weighted normalization, obtain the average confidence level of the line. For each line, take the average confidence level Conf = 1–λ_min / λ_max within a 1 m buffer zone. A value greater than 0.7 is considered high confidence, between 0.4 and 0.7 is considered medium confidence, and less than 0.4 is considered low confidence. Step S52: Perform least-squares straight-line fitting on the first and last nodes of the trailing edge in 3D space, project it onto the horizontal plane, calculate the clockwise angle with due north, and automatically output the 0-360° direction value; the 3D coordinates of the first and last nodes are (x1,y1,z1)→(x2,y2,z2), and the angle between the horizontal projection and due north is Azi = arctan2(Δx,Δy)×180 / π, which is used to calculate the direction; Step S53: Determine the position of the trailing edge based on topological relationships. If it is located at the outer edge of the top of the landslide body, it is marked as the outer trailing edge; if it is located at the inner edge, it is marked as the inner trailing edge; if it is isolated on one side, it is marked as pending verification. Complete the assignment of geometric type attributes; make topological distinction with the 1 m buffer zone of the adjacent landslide surface to determine the geometric type.

[0014] As a preferred option, step S6 specifically includes the following steps: Step S611: For any registration period T k trailing edge line L k Densify the vertex count until the average distance between the vertex and L0 is ≤2.5m; Step S612: For L k Each vertex v k(j) Find the closest point v in L0 by Euclidean distance. 0(i) If the distance is greater than d_max, then discard. Step S613: Construct a joint cost function, taking into account both Euclidean distance minimization and path consistency constraints; Step S621: Register the two adjacent L periods k ′ and L k ₊1′, establish a bidirectional nearest point mapping to ensure a one-to-one correspondence; Step S622: In L k ′ Each point p k (i) A local Frenet frame (tangential t, normal n, and subnormal b) is established at point (i), and the corresponding point p is obtained by interpolation along the normal n at a grid resolution of 2.5 m. k ₊1(i); Step S623: Calculate the three-dimensional displacement vector, the change in orientation, and the change in curvature; Step S624: Synchronously record the confidence weight of this point; Step S625: Output the difference feature set F k ={Δd, Δθ, Δκ, w}, and store them in the spatiotemporal difference database.

[0015] As a preferred option, step S7 specifically includes the following steps: Step S711: The time window length is fixed at n=5 periods; Step S712: For each point i, take the set of neighboring points N(i) along the arc length of the trailing edge ±50 m for the spatial window length; Step S713: Weight Model W t = λ^(T_target – T k ), λ=0.8, W s = exp(–‖s–s_i‖² / 2σ²), σ=25 m, W_c = w(i), and the comprehensive weight W(i,k)=W t ·W s ·W_c / ΣW; Step S714: For target point i during the prediction period T₊1, the increment Δd̂(i) = ΣkΣ_{j∈N(i)} W(j,k)·Δd k (j), in the same way, Δθ̂(i) and Δκ̂(i) are obtained, forming the prediction feature F̂={Δd̂, Δθ̂, Δκ̂}; Step S721: Shift each point p_last(i) of the latest L_last by Δd̂(i) to obtain the coarse prediction point p̂(i); Step S722: Use cubic B-spline smoothing, control point spacing ≤ 10 m, curvature constraint |κ| ≤ 1 / 50 m⁻¹, to eliminate extrapolation noise; Step S723: Perform 3D Kuwahara filtering on the smoothed prediction line again to retain significant edges and obtain the final predicted trailing edge line L_pred.

[0016] By employing the above technical solutions, this invention has the following beneficial effects compared to existing technologies: 1. The entire process in a single phase does not require a slope threshold, overcoming the bottlenecks of inaccurate calculations and difficult parameter tuning in traditional methods; 2. Relying solely on multibeam DEMs eliminates the need for additional seismic profiles, significantly reducing survey costs; 3. The entire process eliminates the need for manual delineation of potential slip zones, solving the interpretability problem; 4. The algorithm has strong noise resistance and is stable at resolutions of 2.5 m, 5 m, and 7.5 m, making it suitable for different water depths in nearshore waters worldwide; 5. The process from single-stage processing to predictive output can be completed in a short time, meeting the real-time monitoring needs of marine engineering.

[0017] Additional aspects and advantages of the invention will become apparent in the following description or may be learned by practice of the invention. Attached Figure Description

[0018] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which: Figure 1 The overall flowchart of intelligent identification and deformation prediction of submarine landslides based on time-series multibeam bathymetry data is as follows: It shows the process of intelligent extraction of the trailing edge line for a single period by acquiring raw multibeam bathymetry data, sequentially performing GST coherence weighting, 3D-WPCA, Kuwahara filtering and region growing algorithm; then realizing multi-period vector space registration by improving ICP, and automatically generating the predicted trailing edge line for the next period by using spatiotemporal sliding window TSW weighted extrapolation. Detailed Implementation

[0019] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0020] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and therefore the scope of protection of the invention is not limited to the specific embodiments disclosed below.

[0021] The following is combined Figure 1 The present invention provides a detailed description of the intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to embodiments of the present invention.

[0022] like Figure 1 As shown, this invention provides a method for intelligent identification and deformation prediction of submarine landslides based on single-period intelligent extraction, multi-period spatial registration, and spatiotemporal model prediction. The method is characterized by the following steps: Step S1: Acquisition and preprocessing of raw data in multiple periods T0, T1, and T2. The target sea area is measured using the same shipborne multibeam system to obtain raw data. After acoustic ray tracking, tide level correction, and noise reduction, a 2.5 m × 2.5 m regular grid DEM is output, namely DEM_T0, DEM_T1, and DEM_T2, with an elevation accuracy better than 0.3 m. Step S2: Calculate the coherence of the DEM for each single period using the Gradient Structure Tensor (GST) to obtain the reliability weights of each grid point; the GST coherence calculation specifically includes the following steps: Step S21: Calculate ∂h / ∂x, ∂h / ∂y, and ∂h / ∂z pixel by pixel for DEM_Ti to construct the gradient vector ∇h; Step S22: Accumulate the outer product in the 3×3×3 neighborhood to obtain the gradient structure tensor: GST = Σ (∇h)(∇h)^T; Step S23: Find the eigenvalues ​​λ1≥λ2≥λ3, and calculate coherence_GST = (λ1−λ3) / λ1, taking values ​​from 0 to 1. The closer to 1, the more uneven the local surface is, which is used as the reliability weight for subsequent 3D-WPCA.

[0023] Step S3: Using the coherence from Step S2 as weights, perform 3D weighted principal component analysis (3D-WPCA) within a 5×5×5 window to obtain the gradient and dip angle, generating an initial landslide trailing edge probability volume. Apply 3D Kuwahara filtering to the probability volume, preserving edges and suppressing noise based on the principle of minimizing the variance of (2n+1)×(2n+1) sub-windows, to obtain an enhanced landslide trailing edge probability map. Automatically select seed points on the probability map and execute a region growing algorithm to aggregate similar probability pixels into continuous linear objects, completing the initial extraction of the landslide trailing edge. Specifically, this includes the following steps: Three-dimensional weighted principal component analysis includes the following steps: Step S311: Run WPCA point-by-point on DEM_Ti using a 5×5×5 voxel sliding window (actual area ≈ 12.5 m × 12.5 m × 5 m), where C = Σ w · (∇h)(∇h)^T, weight w = exp(−d / d0) · coherence_GST, d0=2 voxels; Step S312: Extract the minimum eigenvalue λ_min of C, which is defined as the roughness index U; U greater than 0.65 is considered as a potential trailing edge voxel; Step S313: The output is a three-dimensional probability volume Probability_Ti, with a value of 0-1. The larger the value, the more likely it is to be the trailing edge of the landslide. The three-dimensional Kuwahara filter preserves edges and suppresses noise, specifically including the following steps: Step S321: Divide Probability_Ti into 8 overlapping 3×3×3 sub-blocks; Step S322: Calculate the mean m_k and standard deviation σ_k of each sub-block, and retain the mean of the sub-block corresponding to σ_min as the filter output; Step S323: Preserve the true trailing edge abrupt change signal and suppress side noise.

[0024] Region-based clustering, specifically including the following steps: Step S331: Automatic seed selection: Find connected cores with local maxima and U>0.75 within the filtered U-body; Step S332: Use the similarity criterion |U_i – U_seed|<0.1 to determine whether to retain a continuous number of voxels ≥200 m² (≈32 voxels); Step S333: The output is a raster cluster Cluster_Ti, with each cluster assigned a unique value.

[0025] Step S3 requires no manual setting of the slope threshold and maintains recognition stability for DEMs with resolutions of 2.5 m, 5 m, and 7.5 m.

[0026] Step S4: Reduce the dimensionality of the 3D probabilistic volume to a 2D plane, and use a feature enhancement algorithm to normalize and enhance the probabilistic map, causing the pixel values ​​at the landslide trailing edge to become polarized. After skeletonization, broken line connection, and spur removal, the enhanced binary image is processed to obtain a single-pixel wide continuous landslide trailing edge line. Vectorization tracking is performed along the skeleton line to generate independent vector line segments, which are then reprojected onto the original 3D DEM space according to spatial coordinates, realizing the vectorization extraction and 3D mapping of each submarine landslide trailing edge. Specifically, this includes the following steps: Step S41: Project Cluster_Ti onto a two-dimensional plane with the highest probability along the depth direction; Step S42: Use grayscale remapping I′ = (I–I_min)^γ / (I_max–I_min)^γ, γ=0.5 to achieve pixel polarization; Step S43: First, binarize the image to a threshold of 0.5 to obtain a binary image B. After refining image B into a single-pixel wide skeleton, automatically connect the broken ends with a Hausdorff distance of less than 3 pixels and a direction difference of less than 30°. Automatically remove burrs with a length of less than 5 pixels. Only retain single trailing edge lines with a length of 40 m or more.

[0027] Step S5: Automatically calculate the confidence level, direction, and geometric type attributes of each vectorized landslide trailing edge vector to achieve refined extraction of the landslide trailing edge; specifically including the following steps: Step S51: Map the vector of the trailing edge of a single landslide back to the original 3D DEM. Extract the GST coherence, 3D-WPCA gradient magnitude, and Kuwahara variance within the buffer zone. After weighted normalization, obtain the average confidence level of the line. For each line, take the average confidence level Conf = 1–λ_min / λ_max within a 1 m buffer zone. A value greater than 0.7 is considered high confidence, between 0.4 and 0.7 is considered medium confidence, and less than 0.4 is considered low confidence. Step S52: Perform least-squares straight-line fitting on the first and last nodes of the trailing edge in 3D space, project it onto the horizontal plane, calculate the clockwise angle with due north, and automatically output the 0-360° direction value; the 3D coordinates of the first and last nodes are (x1,y1,z1)→(x2,y2,z2), and the angle between the horizontal projection and due north is Azi = arctan2(Δx,Δy)×180 / π (0–360°, retain 1 decimal place), which is used to calculate the direction; Step S53: Determine the position of the trailing edge based on topological relationships. If it is located at the outer edge of the top of the landslide body, it is marked as the outer trailing edge; if it is located at the inner edge, it is marked as the inner trailing edge; if it is isolated on one side, it is marked as pending verification. Complete the assignment of geometric type attributes; make topological distinction with the 1 m buffer zone of the adjacent landslide surface to determine the geometric type.

[0028] Step S6: Using the initial landslide trailing edge vector as the primary reference, the improved ICP algorithm is used to complete the spatial registration of multi-stage trailing edge lines, and the displacement, strike, and curvature differences between adjacent stages are calculated point by point along the normal. Specifically, this includes the following steps: Step S61: Using the initial landslide trailing edge vector as the main reference, the improved ICP algorithm is used to spatially register the trailing edge lines of subsequent phases. The registration cost function simultaneously considers the minimization of Euclidean distance and the constraint of strike consistency. Step S611: For any registration period T k The trailing edge L (k≥1) k Densify the vertex count to ensure that the average distance between the vertex and L0 is ≤2.5 m; Step S612: For L k Each vertex v k(j) Find the closest point v in L0 by Euclidean distance. 0(i) If the distance is greater than d_max (initially d_max = 15 m, decreasing by 10% per iteration), then discard the method. Step S613: Construct a joint cost function, taking into account both Euclidean distance minimization and path consistency constraints; Step S62: After registration, calculate the displacement vector between adjacent periods point by point along the normal direction of the trailing edge with a grid resolution of 2.5 m, and simultaneously extract the change in orientation and the change in curvature. Step S621: Register the two adjacent L periods k ′ and L k ₊1′, establish a bidirectional nearest point mapping to ensure a one-to-one correspondence; Step S622: In L k ′ Each point p k (i) A local Frenet frame (tangential t, normal n, and subnormal b) is established at point (i), and the corresponding point p is obtained by interpolation along the normal n at a grid resolution of 2.5 m. k ₊1(i); Step S623: Calculate the three-dimensional displacement vector, the change in orientation, and the change in curvature; Step S624: Synchronously record the confidence weight of this point; Step S625: Output the difference feature set F k ={Δd, Δθ, Δκ, w}, and store them in the spatiotemporal difference database.

[0029] Step S7: Introduce a spatiotemporal sliding window (TSW) to perform a weighted average of the difference features for each single period. The weights are determined by the confidence level and the time distance from the prediction period, resulting in predicted difference features. These predicted difference features are then superimposed onto the latest period's trailing edge, and extrapolated point-by-point to generate the trailing edge for the prediction period. This achieves submarine landslide deformation prediction (i.e., prediction of the spatial displacement, direction, and curvature changes of the landslide trailing edge), completing the intelligent identification and deformation prediction of submarine landslides from time-series multibeam data (multi-period data). Specifically, this includes the following steps: Step S71: Introduce a spatiotemporal sliding window (TSW) to perform a weighted average of the difference features for each single period. The weights are determined by the confidence level and the time distance from the prediction period to obtain the predicted difference features. Step S711: The time window length is fixed at n=5 periods (if the historical period is less than 5 periods, then all periods are taken); Step S712: For each point i, take the set of neighboring points N(i) along the arc length of the trailing edge ±50 m for the spatial window length; Step S713: Weight Model W t = λ^(T_target – T k ), λ=0.8 (time decay), W s = exp(–‖s–s_i‖² / 2σ²), σ=25 m (spatial Gaussian kernel), W_c = w(i) (confidence level), and the comprehensive weight W(i, k )=W t ·W s ·W_c / ΣW; Step S714: The increment Δd̂(i) of target point i during the prediction period T₊1 is Σ k Σ_{j∈N(i)} W(j,k)·Δd k (j), in the same way, Δθ̂(i) and Δκ̂(i) are obtained, forming the prediction feature F̂={Δd̂, Δθ̂, Δκ̂}; Step S72: Overlay the predicted difference features onto the latest trailing edge line, and extrapolate point by point to generate the trailing edge line for the prediction period, so as to realize intelligent prediction of submarine landslide deformation based on multi-period time-series multibeam data. Step S721: Shift each point p_last(i) of the latest L_last by Δd̂(i) to obtain the coarse prediction point p̂(i); Step S722: Use cubic B-spline smoothing, control point spacing ≤ 10 m, curvature constraint |κ| ≤ 1 / 50 m⁻¹, to eliminate extrapolation noise; Step S723: Perform a third-dimensional Kuwahara filter (3×3×3 window) on the smoothed prediction line to preserve significant edges and obtain the final predicted trailing edge line L_pred.

[0030] Terminology Explanation: In this invention, "temporal multibeam data" refers to multi-period temporal multibeam data; "deformation prediction" refers to submarine landslide deformation prediction; "multi-period" refers to multibeam data acquired at two or more time points; "GST" is Gradient Structure Tensor; "3D-WPCA" is Three-Dimensional Weighted Principal Component Analysis, which is used for description only and does not constitute a limitation of the scheme.

[0031] In the description of this specification, the terms "one embodiment," "some embodiments," "specific embodiment," etc., refer to a specific feature, structure, material, or characteristic described in connection with that embodiment or example, which is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0032] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for intelligent identification and deformation prediction of submarine landslides based on time-series multibeam data, characterized in that... Specifically, it includes the following steps: Step S1: Acquisition and preprocessing of raw data in multiple periods T0, T1, and T2. The target sea area is measured using the same shipborne multibeam system to obtain raw data. After completing acoustic ray tracking, tide level correction, and noise reduction, a 2.5 m × 2.5 m regular grid DEM is output, namely DEM_T0, DEM_T1, and DEM_T2. Step S2: Calculate the coherence of the DEM for each single period using the gradient structure tensor GST to obtain the reliability weight of each grid point. Step S3: Using the coherence from Step S2 as weights, perform 3D weighted principal component analysis (3D-WPCA) within a 5×5×5 window to obtain the gradient and dip angle, generating an initial landslide trailing edge probability volume; apply 3D Kuwahara filtering to the probability volume, preserving edges and suppressing noise according to the principle of minimizing the variance of (2n+1)×(2n+1) sub-windows, to obtain an enhanced landslide trailing edge probability map; automatically select seed points on the probability map, execute a region growing algorithm, and aggregate similar probability pixels into continuous linear objects to complete the initial extraction of the landslide trailing edge; Step S4: Reduce the dimensionality of the 3D probabilistic volume to a 2D plane, and use a feature enhancement algorithm to normalize and enhance the probabilistic map, making the pixel values ​​at the landslide trailing edge polarized; after skeletonization, broken line connection and spur removal, the enhanced binary image is processed to obtain a single-pixel wide continuous landslide trailing edge line; vectorization tracking is performed along the skeleton line to generate independent vector line segments, and the segments are reprojected onto the original 3D DEM space according to the spatial coordinates, realizing the vectorization extraction and 3D mapping of each submarine landslide trailing edge; Step S5: Automatically calculate the confidence level, direction and geometric type attributes of each vectorized landslide trailing edge vector to achieve refined extraction of landslide trailing edge; Step S6: Using the trailing edge vector of the first landslide as the main reference, the improved ICP algorithm is used to spatially register the trailing edge lines of subsequent phases. The registration cost function considers both Euclidean distance minimization and strike consistency constraints. After registration, the displacement vector between adjacent phases is calculated point by point along the normal direction of the trailing edge line at a grid resolution of 2.5 m, and the strike change and curvature change are extracted simultaneously. Step S7: Introduce a spatiotemporal sliding window (TSW) to perform a weighted average of the difference features for each single period. The weights are determined by the confidence level and the time distance from the prediction period to obtain the predicted difference features. The predicted difference features are then superimposed on the latest period trailing edge line, and the prediction period trailing edge line is generated by extrapolating point by point to achieve intelligent prediction of multi-period multi-beam submarine landslides.

2. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... The GST coherence calculation in step S2 specifically includes the following steps: Step S21: Calculate ∂h / ∂x, ∂h / ∂y, and ∂h / ∂z pixel by pixel for DEM_Ti to construct the gradient vector ∇h; Step S22: Accumulate the outer product in the 3×3×3 neighborhood to obtain the gradient structure tensor: GST = Σ (∇h)(∇h)^T; Step S23: Find the eigenvalues ​​λ1≥λ2≥λ3, and calculate coherence_GST = (λ1−λ3) / λ1, taking values ​​from 0 to 1. The closer to 1, the more uneven the local surface is, which is used as the reliability weight for subsequent 3D-WPCA.

3. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... The three-dimensional weighted principal component analysis in step S3 specifically includes the following steps: Step S311: Run WPCA point-by-point on DEM_Ti using a 5×5×5 voxel sliding window (actual area ≈ 12.5 m × 12.5 m × 5 m), where C = Σ w · (∇h)(∇h)^T, weight w = exp(−d / d0) · coherence_GST, d0=2voxel; Step S312: Extract the minimum eigenvalue λ_min of C, which is defined as the roughness index U; U greater than 0.65 is considered as a potential trailing edge voxel; Step S313: The output is a three-dimensional probability volume Probability_Ti, with a value of 0-1. The larger the value, the more likely it is to be the trailing edge of the landslide.

4. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 3, characterized in that... The three-dimensional Kuwahara filtering in step S3, which preserves edges and suppresses noise, specifically includes the following steps: Step S321: Divide Probability_Ti into 8 overlapping 3×3×3 sub-blocks; Step S322: Calculate the mean m_k and standard deviation σ_k of each sub-block, and retain the mean of the sub-block corresponding to σ_min as the filter output; Step S323: Preserve the true trailing edge abrupt change signal and suppress side noise.

5. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 4, characterized in that... The region growing clustering in step S3 specifically includes the following steps: Step S331: Automatic seed selection: Find connected cores with local maxima and U>0.75 within the filtered U-body; Step S332: Use the similarity criterion |U_i – U_seed|<0.1 to make a judgment; only those with a continuous number of voxels ≥200 m² are retained. Step S333: The output is a raster cluster Cluster_Ti, with each cluster assigned a unique value.

6. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... Step S4 specifically includes the following steps: Step S41: Project Cluster_Ti onto a two-dimensional plane with the highest probability along the depth direction; Step S42: Use grayscale remapping I′ = (I–I_min)^γ / (I_max–I_min)^γ, γ=0.5 to achieve pixel polarization; Step S43: First, binarize the image to a threshold of 0.5 to obtain a binary image B. After refining image B into a single-pixel wide skeleton, automatically connect the broken ends with a Hausdorff distance of less than 3 pixels and a direction difference of less than 30°. Automatically remove burrs with a length of less than 5 pixels. Only retain single trailing edge lines with a length of 40 m or more.

7. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... Step S5 specifically includes the following steps: Step S51: Map the vector of the trailing edge of a single landslide back to the original 3D DEM. Extract the GST coherence, 3D-WPCA gradient magnitude, and Kuwahara variance within the buffer zone. After weighted normalization, obtain the average confidence level of the line. For each line, take the average confidence level Conf = 1–λ_min / λ_max within a 1 m buffer zone. A value greater than 0.7 is considered high confidence, between 0.4 and 0.7 is considered medium confidence, and less than 0.4 is considered low confidence. Step S52: Perform least-squares straight-line fitting on the first and last nodes of the trailing edge in 3D space, project it onto the horizontal plane, calculate the clockwise angle with due north, and automatically output the 0-360° direction value; the 3D coordinates of the first and last nodes are (x1,y1,z1)→(x2,y2,z2), and the angle between the horizontal projection and due north is Azi = arctan2(Δx,Δy)×180 / π, which is used to calculate the direction; Step S53: Determine the position of the trailing edge based on topological relationships. If it is located at the outer edge of the top of the landslide body, it is marked as the outer trailing edge; if it is located at the inner edge, it is marked as the inner trailing edge; if it is isolated on one side, it is marked as pending verification. Complete the assignment of geometric type attributes; make topological distinction with the 1 m buffer zone of the adjacent landslide surface to determine the geometric type.

8. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... Step S6 specifically includes the following steps: Step S611: For any registration period T k trailing edge line L k Densify the vertex count to ensure that the average distance between the vertex and L0 is ≤2.5 m; Step S612: For L k Each vertex v k(j) Find the closest point v in L0 by Euclidean distance. 0(i) If the distance is greater than d_max, then discard. Step S613: Construct a joint cost function, taking into account both Euclidean distance minimization and path consistency constraints; Step S621: Register the two adjacent L periods k ′ and L k ₊1′, establish a bidirectional nearest point mapping to ensure a one-to-one correspondence; Step S622: In L k ′ Each point p k (i) A local Frenet frame (tangential t, normal n, and subnormal b) is established at point (i), and the corresponding point p is obtained by interpolation along the normal n at a grid resolution of 2.5 m. k ₊1(i); Step S623: Calculate the three-dimensional displacement vector, the change in orientation, and the change in curvature; Step S624: Synchronously record the confidence weight of this point; Step S625: Output the difference feature set F k ={Δd, Δθ, Δκ, w}, and store them in the spatiotemporal difference database.

9. The intelligent identification and deformation prediction method for submarine landslides based on time-series multibeam data according to claim 1, characterized in that... Step S7 specifically includes the following steps: Step S711: The time window length is fixed at n=5 periods; Step S712: For each point i, take the set of neighboring points N(i) along the arc length of the trailing edge ±50 m for the spatial window length; Step S713: Weight Model W t = λ^(T_target – T k ), λ=0.8, W s = exp(–‖s–s_i‖² / 2σ²), σ=25 m, W_c = w(i), and the comprehensive weight W(i,k)=W t ·W s ·W_c / ΣW; Step S714: The increment Δd̂(i) of target point i during the prediction period T₊1 is Σ k Σ_{j∈N(i)} W(j,k)·Δd k (j), in the same way, Δθ̂(i) and Δκ̂(i) are obtained, forming the prediction feature F̂={Δd̂, Δθ̂, Δκ̂}; Step S721: Shift each point p_last(i) of the latest L_last by Δd̂(i) to obtain the coarse prediction point p̂(i); Step S722: Use cubic B-spline smoothing, control point spacing ≤ 10 m, curvature constraint |κ| ≤ 1 / 50 m⁻¹, to eliminate extrapolation noise; Step S723: Perform 3D Kuwahara filtering on the smoothed prediction line again to retain significant edges and obtain the final predicted trailing edge line L_pred.

Images