Mine depth prediction method based on CEL method and multi-year sedimentary evolution data

Abstract

This invention discloses a method for predicting the burial depth of sea mines based on the CEL method and multi-year sedimentary evolution data, belonging to the field of sea mine detection and marine engineering. This method solves the problem that existing technologies struggle to systematically predict the final burial depth of historically abandoned sea mines over large sea areas across time scales. The technical solution involves: establishing a three-dimensional finite element model of the sea mine penetrating the seabed soil; calculating the vertical component of the bottom-touching velocity at different water depths by setting soil constitutive parameters and combining the physical characteristics of the sea mine with entry conditions; and simulating the initial penetration depth of the sea mine. Simultaneously, acquiring multi-period historical nautical chart water depth data of the predicted sea area; and obtaining the long-term erosion and deposition thickness sequence {Δ} through a unified reference surface and comparative analysis. H t The final burial depth distribution for the termination year is obtained recursively based on this sequence. H T and the distribution of the maximum possible burial depth within a historical period H max ; the initial burial depth and {Δ H t The distribution map of predicted mine burial depth is generated by overlaying and recursively extrapolating the data according to the time series.

Description

Technical Field

[0001] This invention relates to the fields of mine detection and marine engineering. More specifically, this invention relates to a method for predicting mine burial depth based on the CEL method and multi-year sedimentary evolution data. Background Technology

[0002] In the field of underwater unexploded ordnance (such as mines) detection and disposal engineering, accurately predicting their burial depth in underwater sedimentary layers is a core prerequisite for scientifically selecting detection technologies, optimizing mine clearance processes, and controlling engineering risks and costs. This need is particularly urgent in the construction of major marine engineering projects such as offshore wind farms in areas with historical minefields. However, existing technological systems still face a series of interconnected and fundamental challenges in achieving effective and reliable prediction of the final burial depth of historical mines over large areas of sea.

[0003] First, from the perspective of overall prediction methodologies, existing technologies often exhibit limitations, failing to provide comprehensive analysis across time scales. On one hand, researchers have developed impact penetration analysis methods, including one-dimensional, two-dimensional, and even three-dimensional models, to calculate the "initial penetration depth" of mines at the moment of impact. These studies primarily focus on transient factors such as mine shape, entry velocity, and seabed soil composition, with time scales typically on the order of seconds. On the other hand, marine geology has long been dedicated to studying the topographic evolution of coastal zones and seabeds. By analyzing historical nautical charts and hydrological sediment data from multiple years, it is possible to reveal the "long-term erosion and deposition patterns" of marine areas over decades. However, these two technical approaches, focusing on "short-term impact" and "long-term evolution" respectively, have long been disconnected. For historical mines deployed decades ago, their final state is the product of the combined effects of the initial impact and the subsequent long-term sedimentary environment. Existing technologies lack a predictive framework that systematically couples these two physical processes with different time dimensions. This makes it impossible to make a macroscopic and reasonable prediction of the potential and maximum burial depth of mines from the engineering planning stage. Consequently, the design of subsequent detection schemes lacks sufficient scientific basis, which may result in either underestimation and missed detection or excessive conservatism and waste of resources.

[0004] Secondly, within the two analytical dimensions mentioned above, existing technologies also have specific bottlenecks affecting prediction accuracy. In the initial penetration depth prediction stage, although numerical simulation methods (such as the finite element method) are widely used, the high-speed penetration of seabed soil by mines is a dynamic challenge involving large deformations, material nonlinearity, and complex contacts. Constructing a simulation model that can maintain computational stability while accurately capturing the physical processes is inherently challenging, including but not limited to the appropriate selection and parameter determination of the soil constitutive model, the handling of fluid-structure interaction boundaries, and the balance between computational scale and efficiency. Furthermore, the calculation of the mine's "bottoming velocity," a key input condition for simulation, is often based on the simplified assumption of vertical entry into the water. However, in actual airdrop operations, the mine may have a certain entry angle, and its horizontal velocity component will affect the penetration attitude and depth. Ignoring this factor may lead to deviations in the assessment of the most unfavorable condition (maximum penetration depth).

[0005] In long-term seabed evolution analysis, current practices largely rely on comparing water depth and topographic data from different periods to obtain net results of scouring and deposition (such as maximum deposition thickness). While this method can reflect the overall trend, it compresses decades of evolution into a static "final state," losing dynamic sequence information about the multiple scouring and backfilling processes that may occur during the process. This "static snapshot" analysis cannot support the tracing or prediction of the complex evolutionary paths that mine burial processes may undergo over time, limiting the depth and flexibility of the analysis.

[0006] Furthermore, existing methods, whether analyzing initial penetration or long-term evolution, typically treat the seabed as a homogeneous system controlled solely by natural hydrodynamics (such as tides and waves). However, in the actual marine environment, submarine engineering activities (such as channel dredging and pipeline laying) and benthic community activities significantly and unconventionally disturb the stability and transport patterns of local sediments. These local erosion and deposition anomalies caused by such "non-natural" driving forces are error sources that traditional prediction models based on hydrodynamic and topographic analysis struggle to identify and quantify. Ignoring these factors can lead to significant deviations between predicted results and actual conditions in specific local areas, but current technologies lack mature methods for effectively integrating such diverse environmental disturbance information.

[0007] Finally, examining the entire chain of engineering applications, the reliability and persuasiveness of a prediction method ultimately require verification through field experiments. Currently, predictions of mine burial depth over large-scale sea areas mostly remain at the level of theoretical calculations or case studies, lacking systematic, large-scale engineering practice-integrated field detection data to quantitatively verify and provide feedback loops for the prediction results. This lack or weakness in the verification process means that the credibility assessment and subsequent optimization of the prediction model lack a solid data foundation, and also affects the degree to which this technology is directly adopted in major engineering decisions.

[0008] In summary, the main challenges currently facing the field of mine burial depth prediction are: macroscopically, a lack of prediction methodologies that couple transient dynamics with long-term geological evolution; microscopically, numerous technical bottlenecks exist in areas such as accurate modeling, key parameter calculation, dynamic process reconstruction, and consideration of complex environmental factors; and the entire prediction process lacks an empirical verification stage closely integrated with engineering practice. These problems make efficient, accurate, and reliable pre-assessment of mine burial conditions in large-scale historical minefields a challenging technical problem. Summary of the Invention

[0009] This invention provides a method for predicting the burial depth of sea mines based on the CEL method and multi-year sedimentary evolution data, comprising the following steps: S1. Establish a three-dimensional finite element model of a mine penetrating the seabed soil. The model is constructed using the Coupled Eulerian-Lagrange (CEL) method, where the mine entity is simulated using a Lagrange mesh and the seabed soil is simulated using an Eulerian mesh. S2. Set the material constitutive model and corresponding parameters of the seabed soil in the three-dimensional finite element model. The parameters include the mass of the mine, the cross-sectional area of ​​the mine, the drag coefficient of the mine, the density of the seabed soil, the internal friction angle, the cohesion, and the friction coefficient between the mine and the soil. S3. Based on the mine's mass, dimensions, drag coefficient, seawater density, and water depth, and considering the mine's initial entry conditions, calculate the vertical component of the mine's bottom-touching velocity corresponding to different water depths. Input this vertical component of the bottom-touching velocity as initial conditions into the three-dimensional finite element model for explicit dynamic analysis, simulating the mine's centroid penetration depth in the seabed soil. H p ; S4. Obtain historical nautical chart depth data for the predicted sea area over a historical period, unify the historical nautical chart depth data to the same depth datum, and obtain the seabed sedimentation thickness sequence {Δ} for the predicted sea area over a historical period through comparative analysis. H t}; S5, Based on the centroid penetration depth H p Determine the initial burial depth of the mine's top based on its external dimensions. H i and will H i With {Δ H t The final burial depth distribution is obtained by recursively stacking the time series data. H T or maximum possible burial depth distribution H max Output a distribution map of the predicted burial depth.

[0010] Preferably, in step S1, the process of establishing the three-dimensional finite element model of the mine penetrating the seabed soil specifically includes: Based on the actual external dimensions of the target mine type, including its total length L With the maximum diameter D A three-dimensional geometric model of the mine entity is established in computer-aided design software; The three-dimensional geometric model of the mine entity is imported into finite element analysis software. In the finite element analysis software, based on the mine's mass M, material properties are assigned to the three-dimensional geometric model of the mine entity, including material density. ρ mine Through formula ρ mine = M / V model The calculation yielded, where V model Let V be the volume of the three-dimensional geometric model; In the finite element analysis software, an Eulerian mesh of the seabed soil is created around the three-dimensional geometric model of the mine entity. The dimensions of the Eulerian mesh in the longitudinal and transverse directions are not less than the maximum diameter of the mine. D 80 times the length of the mine, and its vertical dimension is not less than the total length of the mine. L 40 times; In the finite element analysis software, the contact properties between the mine entity and the seabed soil are defined, specifically by creating a dynamic universal contact pair and assigning the friction coefficient between the mine and the soil to the tangential behavior property of the dynamic universal contact pair. In the finite element analysis software, boundary conditions are applied to the model of the seabed soil, specifically constraining the normal displacement of all sides of the model and constraining the bottom surface of the model. X Directional displacement, Y Directional displacement and Z Directional displacement.

[0011] Preferably, in step S2, the material constitutive model of the seabed soil adopts the Mohr-Coulomb elastoplastic constitutive model; The parameters of the Mohr-Coulomb elastoplastic constitutive model are defined, including the density, internal friction angle, and cohesion of the seabed soil, with the density ranging from 1500 kg / m³. 3 Up to 2100 kg / m 3 The internal friction angle ranges from 25° to 40°, and the cohesion ranges from 5 kPa to 30 kPa.

[0012] Preferably, in step S3, the vertical component of the mine's bottom contact velocity is calculated. V(X) The process is as follows: After a mine enters the water, it experiences forces along its trajectory, including gravity, buoyancy, and fluid resistance proportional to the square of its velocity. Let the mass of the mine be... M The projected area facing the wind is S The drag coefficient is C y Seawater density is ρ w The drainage volume is V dis The acceleration due to gravity is g First, calculate the mine's maximum velocity (terminal velocity) in water. v np : ; in, Mg - ρ w gV dis The net sinking force is caused by the difference between gravity and buoyancy. The net acceleration constant 'a' caused by the difference between gravity and buoyancy is defined as follows, neglecting drag: ; Let the instantaneous velocity of the mine upon entering the water be... v 0, water entry angle is θ ( θ=0 (Indicates vertical entry into water), water depth is X According to the one-dimensional velocity decay law along the trajectory direction, the velocity of a mine along its trajectory at depth is... u(X) satisfy: ; Projecting the velocity along the trajectory onto the vertical direction yields the vertical component of the mine's bottom impact velocity: V(X) = u(X) cos θ.

[0013] The vertical component of the bottoming velocity V(X) As the initial input conditions for the three-dimensional finite element model, explicit dynamic analysis is performed to obtain the centroid penetration depth of the mine into the seabed soil. H p .

[0014] Preferably, in step S4, while acquiring the historical nautical chart depth data, marine benthic organism survey data of the predicted sea area are also acquired; based on the biomass in the survey data, a marine biological activity intensity index is quantified and generated. B ; The marine biological activity intensity index BIt is calculated using the following formula: ; in, B α is the biological activity intensity index; α is the standardized coefficient. W i For the first i The perturbation weight per unit biomass of benthic organisms is determined based on their burrowing ability and activity frequency; A i For the first i Biomass of benthic organisms per unit area.

[0015] Preferably, in step S4, records of seabed engineering activities in the predicted sea area are also obtained; based on the marine biological activity intensity index... B Is it greater than the preset threshold? B t Or whether there are records of submarine engineering construction, in the scouring and silting thickness sequence {Δ H t The corresponding area is marked as a disturbed area, or the disturbed area is marked prominently in the final burial depth prediction distribution map / maximum possible burial depth distribution map.

[0016] Preferably, in step S4, the historical time period is divided into: T A series of consecutive sub-time periods ( T >2), and obtain T +1 historical nautical chart depth data to construct a thickness sequence of scour and sedimentation changes.

[0017] Specifically: the historical time period is divided into chronological order. T For each adjacent sub-time period, historical nautical chart depth data for the endpoints of each sub-time period is obtained, resulting in a total of [data missing]. T +1 phase water depth dataset{ D 0, D 1,…, D T The water depth data from each period are unified to the same vertical reference plane, and after spatial registration and interpolation to a uniform grid resolution, the thickness of the scouring and deposition change between adjacent periods is calculated: Δ H t = f ( D t-1 , D t ), t =1,2,…, T ; where Δ H t Indicates the first t A grid of thickness variations corresponding to seabed topography changes within a sub-period; specifying Δ H t>0 indicates siltation (seabed rise), Δ H t <0 indicates scouring (seabed lowering). This yields the set of scouring and deposition thickness changes arranged in a time series {Δ H 1,Δ H 2,…,Δ H T}; The set of thicknesses resulting from scouring and silting changes is used for the subsequent step S5 to calculate the burial depth and extract the maximum possible burial depth.

[0018] Preferably, in step S5, the initial penetration depth H p Let be the vertical penetration depth of the mine's center of mass relative to the initial seabed surface (positive downwards), and from this, determine the initial burial depth of the mine's top relative to the initial seabed surface. H i Specifically, we introduce geometric transformation distance. h ,in h The axial distance from the mine's center of mass to its top is preferably determined by measurement using the mine's three-dimensional geometric model; when measurement is not possible, h Desirable h = kL ,in L Let be the total length of the mine, and k The value ranges from 0.45 to 0.65. Therefore, the initial burial depth of the top of the mine is: H i = H p - h .

[0019] Furthermore, the superposition calculation is performed dynamically using a time-series recursive method: based on the initial burial depth. H i As the initial burial depth H 0 The set of scouring and silting change thicknesses arranged according to the aforementioned time series {Δ H 1,Δ H 2,…,Δ H T}, through the recursive formula H t = H t-1 +Δ H t , ( t =1,2,…, T Calculate the burial depth at the end of each sub-period sequentially. H t This allows us to obtain the final burial depth distribution for the termination year. H T ; and according toH max =max{ H 0 , H 1 ,…, H T The maximum possible burial depth distribution within the historical time period is obtained to dynamically simulate the burial depth evolution process of mines at any point in time within the historical time period and output the final result.

[0020] Preferably, it also includes: S6. Conduct on-site detection and locate N suspected mine targets within the predicted sea area, where N is an integer greater than 10; excavate K targets among the N suspected mine targets, where K is an integer greater than 1 and less than or equal to N; measure the actual burial depth data of the K targets, including the measured burial depth value of each target; compare and analyze the measured burial depth values ​​of the K targets with the predicted burial depth values ​​at the corresponding coordinate positions in the final burial depth prediction distribution map.

[0021] The present invention offers at least the following beneficial effects: The mine burial depth prediction method based on the CEL method and multi-year sedimentary evolution data provides a comprehensive prediction framework that combines the impact dynamics of a mine at the moment of deployment with the natural sedimentary evolution of the seabed over the following decades. This allows the assessment of the final state of historically abandoned mines to take into account both the instantaneous impact effect and the long-term depositional impact, thus overcoming the limitations of previous methods that only focused on a single time scale. It provides a more comprehensive and reliable macro-level decision-making basis for formulating scientific and reasonable overall detection and disposal plans in the early stages of large-scale mine clearance projects.

[0022] By meticulously specifying every step from the three-dimensional geometric reconstruction of the mine to the establishment of the complete finite element model, including determining reasonable model dimensions, assigning accurate material properties, defining precise contact relationships, and setting stable boundary conditions, the constructed numerical model was ensured to effectively and stably simulate the complex and large deformation process of a mine penetrating the soil. This laid a solid foundation for obtaining reliable initial penetration depth calculation results.

[0023] A specific elastoplastic constitutive model, widely validated in geotechnical engineering, is explicitly adopted to characterize seabed soil, and reasonable ranges for its key mechanical parameters are set based on actual engineering geological conditions. This makes the mechanical response of the soil in the simulation closer to reality, significantly improving the engineering credibility and accuracy of the numerical prediction results of the initial penetration depth.

[0024] When calculating the crucial initial condition of the mine's bottom-touching velocity, this study not only considered the mine's basic physical properties and water depth but also innovatively introduced the entry angle parameter. The initial velocity was decomposed into a vertical component using a physical formula. This improvement allows the calculation to more realistically reflect the actual possible inclined entry conditions, thus providing more accurate and diverse input conditions for subsequent penetration simulations and facilitating the assessment of penetration depths in a wider range of scenarios.

[0025] In analyzing historical seabed evolution, in addition to traditional water depth and topographic data, two important local disturbance factors—marine biological activity and seabed human engineering activities—are also considered. By setting thresholds or recording judgments to identify specific areas affected by these factors, the prediction model can identify anomalies beyond natural hydrodynamic laws, thereby improving the overall adaptability and local accuracy of the prediction results in complex real-world environments.

[0026] The aforementioned concept of biological activity intensity is presented with a concrete and feasible quantitative implementation path, namely, by conducting standardized benthic organism surveys and calculating based on the survey data according to predetermined rules. This transforms the improved feature from a theoretical concept into repeatable and operable technical steps, enhancing the practicality and feasibility of the entire method.

[0027] Based on the identification of disturbed areas, a post-processing mechanism is further proposed to downgrade the credibility of prediction results for relevant areas or add uncertainty explanations. This is equivalent to providing clear "accuracy annotations" for the prediction map, which can effectively guide engineers to focus on areas with high prediction uncertainty, thereby taking more prudent or supplementary measures when formulating specific work plans, and enhancing the value of the entire prediction method in engineering risk management.

[0028] By dividing long-term historical periods into multiple consecutive sub-periods and performing detailed erosion and deposition calculations for each adjacent period, a dataset of erosion and deposition changes arranged in a time series is obtained. This method changes the previous approach of only obtaining a final cumulative thickness, revealing the dynamic evolution of the seabed over decades and providing a rich data foundation for more in-depth time-dimensional analysis.

[0029] Based on the aforementioned time-series scouring and sedimentation data, a dynamic simulation of the burial depth of mines after deployment over time was achieved, rather than simply providing a final state. This allows analysts to trace the possible burial depth of mines at any historical moment or analyze their burial rate, greatly expanding the analytical capabilities and application flexibility of this method in historical research, process review, and prediction of future trends.

[0030] A field verification step was added at the end of the prediction process. By conducting actual marine exploration and sampling excavation, and comparing and analyzing the measured data with the prediction results, a complete technical closed loop was formed. This step not only objectively evaluates and verifies the accuracy of the prediction method, but the feedback data generated can also be used for subsequent calibration and optimization of the prediction model, thereby continuously improving the reliability and engineering applicability of the method and forming a virtuous cycle of continuous improvement.

[0031] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0032] Figure 1 Numerical model diagrams of the external structure of the MK13 and MK26 bottom mines.

[0033] Figure 2 This is a three-dimensional model of the mine and the seabed soil.

[0034] Figure 3 This is the mesh generation diagram for the overall finite element model.

[0035] Figure 4 This is a curve showing the relationship between the mine's bottom contact velocity and water depth.

[0036] Figure 5 The graph shows the penetration depth and velocity of MK13 and MK26 mines as a function of time.

[0037] Figure 6 Schematic diagrams of calculation models for vertical penetration (a) and inclined penetration (b) of sea mines.

[0038] Figure 7 This is a curve showing the relationship between the penetration depth and the bottom velocity of a mine.

[0039] Figure 8 This is an underwater topographic map of the sea area studied in 1946.

[0040] Figure 9 This is an underwater topographic map of the sea area studied in 1969.

[0041] Figure 10 This is an underwater topographic map of the sea area studied in 1986.

[0042] Figure 11 This is an underwater topographic map of the sea area studied in 2012.

[0043] Figure 12 This is a diagram showing the changes in seabed erosion and deposition from 1946 to 1969.

[0044] Figure 13This is a diagram showing the changes in seabed erosion and deposition from 1969 to 1986.

[0045] Figure 14 This is a map showing the changes in seabed erosion and deposition from 1986 to 2012.

[0046] Figure 15 This is a map showing the distribution of maximum seabed sediment thickness from 1946 to 2012.

[0047] Figure 16 To study the predicted distribution map of the final burial depth of sea mines (red box indicates areas with burial depth > 4m). Detailed Implementation

[0048] The present invention will now be described in further detail with reference to specific embodiments, so that those skilled in the art can implement it based on the description.

[0049] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0050] It should be noted that, unless otherwise specified, the experimental methods described in the following implementation plan are all conventional methods, and the reagents and materials described are all commercially available unless otherwise specified.

[0051] In the description of this invention, the terms "lateral", "longitudinal", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", and "outer" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0052] This invention provides a method for predicting the burial depth of sea mines based on the CEL method and multi-year sedimentary evolution data, comprising the following steps: S1. Establish a three-dimensional finite element model of a mine penetrating the seabed soil. The model is constructed using the coupled Eulerian-Lagrange method, where the mine entity is simulated using a Lagrange mesh and the seabed soil is simulated using an Eulerian mesh. S2. Set the material constitutive model and corresponding parameters of the seabed soil in the three-dimensional finite element model. The parameters include the mass of the mine, the cross-sectional area of ​​the mine, the drag coefficient of the mine, the density of the seabed soil, the internal friction angle, the cohesion, and the friction coefficient between the mine and the soil. S3. Based on the mine's mass, dimensions, drag coefficient, seawater density, and water depth, and considering the mine's initial entry conditions, calculate the vertical component of the mine's bottom-touching velocity corresponding to different water depths. Input this vertical component of the bottom-touching velocity as initial conditions into the three-dimensional finite element model for explicit dynamic analysis, simulating the initial penetration depth H of the mine in the seabed soil. p ; S4. Obtain the water depth data of multiple historical nautical charts of the predicted sea area within a historical time period, unify the water depth data of the multiple historical nautical charts to the same depth datum, and obtain the cumulative erosion and deposition thickness distribution data of the seabed of the predicted sea area within the historical time period through comparative analysis. S5, Based on the initial penetration depth H p The initial burial depth H is determined by the shape and dimensions of the mine. i The initial burial depth H i The data is overlaid with the cumulative siltation thickness distribution data to generate the final burial depth prediction distribution map of the mines in the predicted sea area.

[0053] In the above technical solution, the modeling and parameter setting section can utilize the widely used commercial finite element analysis software ABAQUS. The three-dimensional solid model of the mine can be created in computer-aided design software, such as using SolidWorks to draw a simplified mine outline. After the model is built, it is imported into ABAQUS using the software's standard data exchange format. In ABAQUS, mine components can be designated as discrete rigid bodies or assigned elastic material properties, with their density calculated based on the mine's total mass and model volume. The soil portion used to simulate the seabed is created using Eulerian mesh technology within the software. The mesh domain's dimensions in both the horizontal and depth directions must be significantly larger than the mine's size, for example, set to tens of times the mine's diameter. The soil's material behavior can be achieved using the Mohr-Coulomb elastoplastic model from the software's material library, with its density, internal friction angle, and cohesion parameters set by referring to typical soil sample test data from the engineering geological survey report of the predicted sea area. The contact properties between the mine surface and the soil are defined as general contact in the software's interaction module, with the friction coefficient given based on material surface property manuals or empirical values. The boundary conditions of the model are set to constrain the normal movement of all sides of the soil region and to completely fix the bottom.

[0054] In the section on calculating the bottoming velocity and simulating penetration, the calculation of the mine's motion parameters after entry into the water can be performed using numerical calculation software such as MATLAB. The required mine mass is calculated. M Geometric dimensions, projected area facing the airflow S drag coefficient C y and drainage volume V disSeawater density can be obtained from historical archives, technical manuals, or measurements taken from a three-dimensional geometric model of a mine. ρ w A standard value can be taken (e.g., 1025 kg / m³). 3 Initial velocity of the mine upon entry into the water. v 0 and water entry angle θ The parameters can be set based on typical airdrop conditions or ballistic simulation results. Based on the above input parameters, the maximum velocity of the mine in water is first calculated. And define the net acceleration constant generated by the difference between gravity and buoyancy when drag is ignored. Subsequently, based on the velocity decay law, the velocity of the mine along its trajectory at a series of representative depth points X within the predicted sea area was calculated. The vertical component of the velocity upon impact with the seabed mud surface is obtained through projection. V(X) = u(X) cos θ The water depth and the vertical component of the velocity { X , V(X) The initial conditions for mine contact with the seabed are input into the finite element model, and explicit dynamic calculations are performed in the ABAQUS / Explicit analysis module. The software dynamically solves the mine's penetration process in the soil based on the set material constitutive, contact, and boundary conditions until it stops. After the calculation, the coordinates of the mine's centroid are extracted using post-processing functions, and its vertical displacement relative to the initial seabed surface is calculated. This vertical displacement is the simulated centroid penetration depth. H p .

[0055] In the seabed evolution analysis and prediction map generation section, it is necessary to collect historical nautical chart depth data for multiple periods in the prediction area. This data can be obtained from maritime or waterway surveying departments. Processing can be done using geographic information system software such as ArcGIS. First, the nautical charts from different periods are georegistered and the depth information is digitized. A key step is to uniformly correct the depth data of all periods to the same vertical datum (e.g., local mean sea level or lowest astronomical tide level), and then perform spatial registration and unified rasterization. The historical time period is divided into... T 1. Continuous sub-time periods, and obtain the corresponding endpoints. T +1 phase water depth dataset{ D 0, D 1,…, D T Subsequently, spatial analysis tools were used to calculate the seabed erosion and deposition thickness Δ between adjacent periods. H t = f ( D t-1 , Dt ), t =1,2,…, T ; forming a set of scouring and silting thickness changes arranged in time series {Δ H 1,Δ H 2,…,Δ H T}, where Δ is specified H t >0 indicates siltation, Δ H t <0 indicates scouring.

[0056] Finally, the superimposed predictions and results are output. First, the centroid penetration depth is obtained from the explicit dynamic analysis. H p The initial burial depth of the mine's top is determined by its external dimensions. H i Introducing geometric transformation distance h ,in h The axial distance from the mine's center of mass to its top is preferably determined by measurement using the mine's three-dimensional geometric model; when measurement is not possible, it can be taken as... h = kL ,in L Let be the total length of the mine, and k The value range is from 0.45 to 0.65. Therefore... H i = H p - h ;by H 0 = H i As the initial burial depth, the thickness is recursively superimposed according to the aforementioned scouring and silting change thickness sequence: H t = H t-1 +Δ H t , ( t =1,2,…, T The final burial depth distribution for the termination year is obtained. H T Meanwhile, to obtain the most probable burial depth distribution over a historical period, the following was extracted based on the maximum value per pixel: H max =max{ H 0 , H 1 ,…, H TThe above overlay, recursion, and maximum value extraction can be completed in batches using the map algebra function of GIS software, ultimately outputting a depth prediction distribution map covering the predicted sea area (including the final depth distribution map). H T / or Maximum possible burial depth distribution map H max To verify the effectiveness of the prediction, the prediction results can be compared and analyzed with the actual mine burial depth data obtained through sonar scanning, magnetic detection, and sampling excavation during the project implementation phase.

[0057] This method, by combining transient dynamic numerical simulation with long-term geological sedimentary data analysis, can provide a relatively comprehensive pre-assessment of the burial status of historical mines over a large area of ​​sea. Its advantages lie in the fact that the use of mature and reliable commercial software and clear historical data as technical support makes the prediction process standardized and repeatable. By coupling the initial impact effect with long-term sedimentation, the prediction results can more reasonably reflect the actual state of the mines after many years, providing a quantitative reference for subsequent detection scheme design, equipment selection, and operational risk assessment, thus helping to improve the scientific rigor and economic efficiency of the early planning of mine clearance projects.

[0058] In other technical solutions, step S1, the process of establishing the three-dimensional finite element model of the mine penetrating the seabed soil, specifically includes: Based on the actual external dimensions of the target mine type, including its total length L With the maximum diameter D A three-dimensional geometric model of the mine entity is established in computer-aided design software; The three-dimensional geometric model of the mine entity is imported into finite element analysis software; in the finite element analysis software, based on the mine's mass... M The three-dimensional geometric model of the mine entity is assigned material properties, including material density. ρ mine Through formula ρ mine = M / V model The calculation yielded, where V model Let V be the volume of the three-dimensional geometric model; In the finite element analysis software, an Eulerian mesh of the seabed soil is created around the three-dimensional geometric model of the mine entity. The dimensions of the Eulerian mesh in the longitudinal and transverse directions are not less than the maximum diameter of the mine. D 80 times the length of the mine, and its vertical dimension is not less than the total length of the mine. L 40 times; In the finite element analysis software, the contact properties between the mine entity and the seabed soil are defined, specifically by creating a dynamic universal contact pair and assigning the friction coefficient between the mine and the soil to the tangential behavior property of the dynamic universal contact pair. In the finite element analysis software, boundary conditions are applied to the model of the seabed soil, specifically constraining the normal displacement of all sides of the model and constraining the bottom surface of the model. X Directional displacement, Y Directional displacement and Z Directional displacement.

[0059] In the above technical solution, establishing a three-dimensional finite element model of a mine penetrating the seabed can be carried out as follows: First, based on the technical drawings or measured data of the target mine (such as the MK13 type), its total length is obtained. L With the maximum diameter D Key external dimensions are then determined. Subsequently, a 3D geometric model of the mine entity can be created using computer-aided design software such as SolidWorks or CATIA, based on these dimensional parameters. During modeling, non-critical features such as threads and small protrusions can be appropriately simplified to improve computational efficiency. The completed 3D model can be saved in a common format such as STEP or IGES.

[0060] Next, import the saved 3D geometric model file of the mine into the preprocessing module of finite element analysis software such as ABAQUS or ANSYS. In the software's material property definition section, material properties need to be assigned to this imported geometry. First, use the software's built-in measurement tools to obtain the volume of the 3D geometric model. V model Then, based on the known total mass of the mines... M Through calculation formula ρ mine = M / V model The material density value that the model should have in the simulation is obtained, and this density value and other necessary mechanical parameters (such as elastic modulus) are input into the software to complete the material definition of the mine component.

[0061] Next, a computational domain for simulating the seabed soil was created in the finite element method software. This domain was discretized using an Eulerian mesh, and its extent needed to completely enclose and be significantly larger than the mine geometry. Specifically, the dimensions of the created Eulerian mesh region in both the length and width directions should be no less than the maximum diameter of the mine. D The length of the mine should be 80 times that of the total length of the mine, and its dimension in the depth direction should not be less than the total length of the mine. L40 times. The mesh type can be the eight-node Eulerian element provided by the software. Then, the mechanical interaction between the mine and the seabed soil needs to be defined. In the interaction module of the software, a dynamic universal contact pair can be created, with the outer surface of the mine geometry set as the master face and the Eulerian region containing the soil set as the slave face. In the properties of this contact pair, the tangential behavior needs to be defined, usually by choosing the penalty function friction formula, and the friction coefficient between the mine and the soil (for example, an empirical value between 0.1 and 0.3 depending on the surface smoothness) is assigned to this property. Finally, the necessary boundary conditions are applied to the seabed soil model to simulate actual constraints: on all the nodes on the sides of the Eulerian mesh, the normal displacement perpendicular to that side is constrained; on all the nodes on the bottom surface of the Eulerian mesh, the normal displacement is simultaneously constrained. X , Y , Z Translational displacements in three directions are used to simulate the supporting effect of a rigid base. After completing all the above steps, an initial finite element model is established that can be used for subsequent explicit dynamic analysis.

[0062] In some other technical solutions, in step S2, the material constitutive model of the seabed soil adopts the Mohr-Coulomb elastoplastic constitutive model. The parameters of the Mohr-Coulomb elastoplastic constitutive model are defined, including the density, internal friction angle, and cohesion of the seabed soil, with the density ranging from 1500 kg / m³. 3 Up to 2100 kg / m 3 The internal friction angle ranges from 25° to 40°, and the cohesion ranges from 5 kPa to 30 kPa.

[0063] In the above technical solution, the constitutive model and parameters for the seabed soil material can be set as follows: In the material module of finite element analysis software such as ABAQUS, select "Mohr-Coulomb Plasticity" from the built-in material model library as the constitutive model for the seabed soil. This model can effectively simulate the elastoplastic behavior and strength characteristics of soil under shear.

[0064] The determination of the parameters required for the model depends on the engineering geological survey data of the predicted sea area. Density parameter ρ The value is usually around 1500 kg / m³ 3 Up to 2100 kg / m 3 The specific value needs to be selected based on the soil type (such as clayey silt, medium-coarse sand) and its physical properties under natural conditions as described in the survey report. (Internal friction angle) φThe value ranges from 25° to 40°. For soils with high sand content, the value is usually close to the upper limit of the range; for silty or clayey soils, the value tends to be closer to the lower limit of the range. Cohesion c The value ranges from 5 kPa to 30 kPa. For non-cohesive sandy soils, this value can be taken at the lower limit of the range or close to zero; for silty clay with a certain degree of cohesion, a higher value can be taken. The specific values ​​of these parameters should be calibrated based on the results of geotechnical tests (such as direct shear tests and triaxial tests). When inputting these parameters into the software, the specific density value determined based on the survey data (such as 1900 kg / m³) must be explicitly entered in the corresponding dialog box. 3 The parameters are: the internal friction angle (e.g., 34°) and the cohesion (e.g., 14 kPa). After inputting the parameters, this material property is assigned to the entire Eulerian grid representing the seabed soil, thus completing the mathematical definition of the soil's mechanical properties.

[0065] In other technical solutions, step S3 involves calculating the vertical component of the mine's bottom-touching velocity. V The specific steps are as follows: Based on the quality of the mine M Frontal projection area S drag coefficient C y Seawater density ρ w Drainage volume V dis and gravitational acceleration g Through formula Calculate the maximum speed of a naval mine in water. v np ;in Mg - ρ w gV dis The net sinking force is the force generated by the difference between gravity and buoyancy. The net acceleration constant generated by the net sinking force is defined. At the angle of entry of the mine into the water θ Under approximately constant conditions, given a vertical water depth Z (That is, the vertical depth of a mine's descent from the water surface to the seabed), which is then converted into the descent distance along its trajectory. X = Z / cos θ And calculate the vertical component of the mine's contact velocity at vertical water depth using the following formula: V(Z) = u(X) cos θ ; where the velocity along the trajectory u(X) We obtain it from the following formula: .

[0066] In the above technical solution, the vertical component of the mine's bottom contact velocity is calculated. V The process can be followed as follows. First, it is necessary to prepare the basic parameters required for the calculation. The mass of the mine... M Frontal projection area S and the drag coefficient C y This information can be obtained from the technical specifications manual or historical archives for this type of mine. Seawater density. ρ w Standard values ​​(e.g., 1025 kg / m³) can be used. 3 Discharge volume of a naval mine. V dis The computer-aided design of 3D models can be used, and the volume measurement function of the software can be used to directly obtain the gravitational acceleration. g Take 9.81 m / s 2 .

[0067] The next step is to perform the first calculation, which is to determine the maximum speed of the mine in water. v np This calculation is based on the physical principle that a mine underwater reaches equilibrium when gravity, buoyancy, and fluid resistance are in balance. Using the prepared parameters, the calculation can be performed using the aforementioned limiting velocity formula. The calculation can be completed in software environments with computational capabilities, such as MATLAB, Python, or Excel.

[0068] Then, the second step of the calculation is performed, which involves solving for different vertical water depths. Z Vertical component of bottoming velocity V(Z) This requires the initial state parameters of the mine upon entry into the water, including the initial entry velocity. v 0 and water entry angle θ . v 0 and θ It can be determined based on historical airdrop records or typical ballistic simulation results. During the calculation, for a series of specific water depth points of interest within the predicted sea area (e.g., from 5m to 40m, in 5m intervals), first... X = Z / cos θ Converted to distance along the trajectory X Then, the velocity along the trajectory is calculated using the velocity decay formula. u(X) Finally, calculate V(Z) = u(X) cos θ Ultimately, these paired { Z , V(Z) The data is compiled into a list and used as input data for defining the initial impact conditions of the mine in subsequent finite element simulations.

[0069] The derivation of the formula for calculating the bottoming speed is as follows: (1) Establish the equation of motion along the trajectory. Let the velocity of the mine along its trajectory after entering the water be... u After entering the water, the forces acting on the object along the direction of motion include: gravity Mg and buoyancy. ρ w gV dis And fluid resistance, which is proportional to the square of the velocity. According to Newton's second law: (1) (2) Define the net acceleration constant and the limiting velocity Neglecting the drag term, the net acceleration constant is: (2) When the speed reaches a stable At that time, the limiting velocity is obtained from equation (1): (3) (3) Rewrite the equations of motion in standard form Substituting equations (2) and (3) into equation (1), we get: (4) (4) Convert the time derivative to the spatial derivative and integrate (expressed as distance along the trajectory). Let represent the sinking distance of the mine along its trajectory, then we have: (5) Substituting equation (4) into equation (5), we get: (6) Taking u(0)=v0 when X=0 as the initial condition, integrating equation (6) yields: (7) thereby (8) (5) Geometric relationship and vertical component of bottoming velocity The angle of entry into the water is θ And when it is approximately constant, the vertical sinking depth Z Distance along the trajectory X satisfy (9) The vertical component of the bottom-touching velocity is the projection of the trajectory velocity onto the vertical direction: (10) In the formula: M : Mass of mines (kg);S Projected area facing the airflow (m²) 2 ); C y Drag coefficient; ρ w Seawater density (kg / m³) 3 ); V dis Discharge volume (m³) 3 g: gravitational acceleration (m / s²) 2 ); v 0: Initial velocity upon entering the water (m / s); θ : Water entry angle ( θ=0 (Indicates vertical entry into water); u(X) Velocity along the trajectory (m / s); v np : Limiting velocity (m / s); a: Net acceleration constant (m / s²) 2 ); X Z: Distance of sinking along the trajectory (m); Z: Vertical water depth (m); V(Z) Vertical component of bottoming velocity (m / s).

[0070] In other technical solutions, in step S4, while acquiring the historical nautical chart depth data, marine benthic organism survey data of the predicted sea area are also acquired; based on the biomass in the survey data, a marine biological activity intensity index is quantified and generated. B ; The marine biological activity intensity index B It is calculated using the following formula: ; in, B α is the biological activity intensity index; W is the standardized coefficient; i The perturbation weight per unit biomass for the i-th type of benthic organism is determined based on its burrowing ability and activity frequency; A i Let be the biomass of the i-th type of benthic organism per unit area.

[0071] In the aforementioned technical solution, the first step in implementation is to obtain marine benthic organism survey data for the predicted sea area. This data can originate from historical ecological survey reports on the area published by maritime authorities, environmental monitoring stations, or research institutions, or can be obtained through supplementary surveys conducted by qualified third-party testing agencies in accordance with standards such as the "Marine Survey Specifications." The survey typically uses box-type or grab-type sediment samples to collect sediment samples from fixed points on the seabed. Through steps such as sieving, sorting, identification, and weighing, a list of benthic organism species and their corresponding biomass data for each sampling station is obtained. The unit of biomass is typically grams per square meter (g / m²). 2The raw data obtained were organized into a structured database containing fields such as station coordinates, survey time, species name, and biomass.

[0072] Next, a marine biological activity intensity index will be generated based on the acquired biomass data. B The specific calculations are performed in data processing software. The required weights are calculated. W i The weighting is determined based on the ecological habits of the target species. For example, for species with deep burrowing depths and frequent activity (such as some bivalves), the weight can be set to 1.2 to 1.5; for species with weaker activity, the weight can be set to 0.5 to 1.0. The specific determination of the weight values ​​can be based on published ecological research literature or by consulting experts in the field. The standardization coefficient α is used to adjust the calculation results to a reasonable numerical range, such as between 0 and 10. Its value can be determined by Σ( ) of all stations. W i · A i The sum of Σ( ) is divided by a preset baseline value to determine the value. In practice, Σ( ) can be calculated from all stations in the surveyed sea area. W i · A i The maximum value of ) is mapped to the exponent. B The upper limit (e.g., 10) is used to calculate the α value. This is done by iterating through the records of each station in the database and substituting them into the formula. B =α·Σ( W i · A i By performing calculations, the index corresponding to each station can be obtained. B .

[0073] Finally, the index calculated from the discrete stations B Spatial interpolation is performed to generate a continuous distribution map covering the entire predicted sea area. This step can be done in Geographic Information System (GIS) software, such as using Kriging interpolation. The software reads data including station coordinates and... B The data table of values ​​is used to generate a raster surface based on spatial correlation rules, namely the marine biological activity intensity index. B A spatial distribution map. This distribution map can serve as a basis for subsequent identification of disturbed areas. For example, a threshold can be set. B t (If based on historical data statistics, B The 70th percentile of the value is set as the threshold, let's say it's 6.5. B Areas with values ​​greater than this threshold are identified as areas of strong biological disturbance.

[0074] This implementation method transforms the qualitative impact of biological activity into a quantitative spatial distribution index through a clear formula and an operable data processing workflow. This provides a repeatable and verifiable technical means for systematically incorporating biological disturbance as an environmental factor into mine burial depth prediction models, helping to improve the adaptability and local accuracy of prediction models in complex real marine environments.

[0075] In some other technical solutions, step S4 also involves acquiring records of seabed engineering activities in the predicted sea area; and based on the marine biological activity intensity index... B Is it greater than the preset threshold? B t Or whether there are records of submarine engineering construction, in the scouring and silting thickness sequence {Δ H t The corresponding area is marked as a disturbed area, or the disturbed area is marked prominently in the final burial depth prediction distribution map / maximum possible burial depth distribution map.

[0076] In the above technical solution, identifying the disturbed area first requires preparing two types of basic data. One is the marine biological activity intensity index calculated as described above. B The spatial distribution of the raster layer is then determined. On the other hand, it is necessary to obtain records of submarine engineering activities in the target historical period from the historical archives of maritime, waterway, or marine engineering construction units. These records are usually in the form of reports, drawings, or databases, and need to be compiled into a list containing information such as project type, construction time, and geographical extent (usually polygon boundary coordinates). Subsequently, in the geographic information system software, the extent of these engineering activities is drawn into a vector polygon layer based on the coordinate information for spatial analysis.

[0077] Next, the disturbed areas are specifically identified. This process is completed in the spatial analysis module of the geographic information system software. First, a threshold for judging biological disturbance is set. B t This threshold can be based on historical data. B Statistical analysis of index data determines, for example, the value to be taken from all historical data. B The 75th percentile of the index value, in this embodiment B t =6.8. Then, using the raster calculator or conditional query tool in the software, for... B The exponential distribution map is used to determine the following conditions: areas with raster cell values ​​greater than 6.8 are initially identified as areas of biological disturbance. Simultaneously, the vector polygon layer representing seabed engineering activities is converted to a format similar to... BThe index map uses raster data with the same spatial reference and resolution, where rasters within the engineering area are assigned a specific identifier value (e.g., 1). Finally, the resulting bio-disturbance identifier raster is logically ORed with the engineering activity identifier raster; that is, if one location satisfies the bio-disturbance strength (…), then… B If any of the conditions in the engineering activities (>6.8) are met, the location is ultimately determined as the "comprehensive disturbed area" and a new binary identifier raster layer is generated, in which the disturbed area is assigned a value of 1 and other areas are assigned a value of 0.

[0078] Ultimately, the resulting "Comprehensive Disturbed Area" marker layer will be used to correct or annotate the seabed scour and deposition analysis results. For example, this marker layer can be overlaid with the "Maximum Sedimentation Thickness Distribution Map," clearly indicating the extent of the disturbed area in the legend. When predicting the final mine burial depth, for the prediction results within the marker area, an explanation can be added to the report, pointing out that due to disturbances by non-natural factors, the scour and deposition patterns in this area are uncertain, and the reliability of the predicted values ​​is relatively low. It is recommended to pay special attention to this area or increase the exploration density in actual engineering surveys. In this way, the prediction results can reflect local anomalies beyond natural sedimentary patterns, improving the overall assessment's detail and the targeted nature of engineering guidance.

[0079] This implementation method establishes a clear and operational spatial identification approach by integrating quantitative biological activity indices with qualitative engineering records. It can systematically identify areas where seabed erosion and deposition patterns deviate from their normal natural state due to intense biological activity or human engineering construction. Separating and treating these areas individually ensures that erosion and deposition analysis and mine burial depth predictions based on historical nautical chart data are no longer uniform, but rather possess spatial differences in "reliability." This provides a more refined and realistic basis for subsequent exploration plan development and engineering risk assessment.

[0080] In some other technical solutions, in step S4, the historical time period is divided into... T A series of consecutive sub-time periods, T The integer is greater than 2; obtain the corresponding value within the historical time period. T +1 historical nautical chart depth data, forming T +1 set of water depth datasets { D 0, D 1,…, D T}; Calculate the seabed erosion and deposition thickness Δ between two adjacent water depth datasets. H t , where is an integer from 1 to , to obtain the set of scouring and silting change thicknesses arranged in time series {Δ H 1,Δ H2,…,Δ H T}

[0081] In the aforementioned technical solution, a more refined temporal analysis of historical seabed evolution can be implemented as follows. First, the total historical period used for analysis (e.g., 1946 to 2012) needs to be divided into consecutive sub-periods. The division can be based on the availability of nautical chart data, key hydrological and sediment events (such as periods of severe storms), or the precision requirements of engineering analysis. For example, based on representative historical nautical chart years, the entire period can be divided into three sub-periods: 1946 to 1969, 1969 to 1986, and 1986 to 2012. T =3.

[0082] Subsequently, historical nautical chart depth data corresponding to the endpoints of each sub-time period were obtained, resulting in a total of T +1 period water depth dataset (in this example, four periods: 1946, 1969, 1986, and 2012, corresponding to { D 0, D 1, D 2, D 3). The water depth data of each period of the nautical chart are scanned (if it is paper), geo-registered, and digital water depth information is extracted in sequence, and uniformly corrected to the same vertical reference plane (such as the local mean sea level or the lowest astronomical tide level). At the same time, spatial registration and uniform rasterization are performed to ensure that the water depth raster of each period has a consistent spatial reference, range and resolution.

[0083] Next, the temporal variation of scour and sedimentation is calculated. Using the spatial analysis tools of Geographic Information System (GIS) software, the differences between adjacent water depth datasets are calculated sequentially over time to obtain the scour and sedimentation variation thickness raster Δ for each sub-period. H t For example, Δ H 1= f ( D 0 , D 1 ) represents the changes from 1946 to 1969, Δ H 2= f ( D 1 , D 2 ) represents the change from 1969 to 1986, Δ H 3= f ( D 2 , D 3 ) represents the changes from 1986 to 2012. Δ is defined as follows. Ht >0 indicates siltation (seabed rise), Δ H t <0 indicates erosion (seabed lowering). This forms a set of erosion and deposition thickness changes arranged in a time series {Δ H 1,Δ H 2,…,Δ H T}, used for subsequent step S5 for recursive calculation of burial depth and extraction of maximum possible burial depth.

[0084] In other technical solutions, in step S5, the superposition calculation is performed dynamically in a recursive manner, specifically: based on the centroid penetration depth obtained from explicit dynamic analysis. H p The initial burial depth of the mine's top relative to the initial seabed surface is determined by the mine's external dimensions. H i and at the initial burial depth H i As the initial burial depth H 0 The set of scouring and silting change thicknesses arranged according to the aforementioned time series {Δ H 1,Δ H 2,…,Δ H T}, through the recursive formula H t = H t-1 +Δ H t , ( t =1,2,…, T Calculate the burial depth distribution at the end of each sub-period sequentially. H t This allows us to obtain the final burial depth distribution for the termination year. H T ; and according to H max =max{ H 0 , H 1 ,…, H T The maximum possible burial depth distribution within a historical time period is obtained, so as to realize the simulation and output of the burial depth evolution of mines at any point in time within the historical time period.

[0085] In the above technical solution, the process of simulating the dynamic evolution of mine burial depth can be operated as follows: Complete the finite element simulation of the initial penetration depth, and the time series set of scour and siltation change thickness {Δ H 1,Δ H 2,…,Δ H TAfter calculation, the relevant raster data are unified into a GIS-compatible format, ensuring that each raster layer has the same spatial extent, projected coordinate system, and cell size. First, the initial raster is obtained through initial burial depth calculation. H 0 And use it as the current burial depth state layer. H current Then, the raster overlay is performed sequentially in chronological order: Step 1: Calculation H 1= H 0 +Δ H 1, and update H current = H 1; Step 2 calculation H 2= H 1 +Δ H 2; and so on, until the result is obtained. H T At the same time, it is possible to { H 0, H 1,…, H T} Perform pixel-by-pixel maximum value extraction to obtain the maximum possible burial depth distribution. H max The above recursive overlay and maximum value extraction can be batch-processed using the map algebra / raster calculator in GIS software, and can be further used to create time series animations or output key year burial depth distribution maps for engineering analysis.

[0086] In some other technical solutions, in step S5, the initial burial depth H i Determined by the following formula: H i = H p - h ;in, H p This represents the vertical penetration depth of the mine's center of mass relative to the initial seabed surface (positive for downward). h The axial distance from the mine's center of mass to its top is preferably determined by measurement using the mine's three-dimensional geometric model; when measurement is not possible, a value is used instead. h = kL ,in L Let be the total length of the mine, and k The value ranges from 0.45 to 0.65.

[0087] In practical implementation, the first step is to obtain the external dimensions of the mine (including its total length). L ( ) and a three-dimensional geometric model. Preferably, the axial distance from the center of mass to the top of the mine is measured directly in the three-dimensional geometric model of the mine. hThis avoids the bias caused by substituting empirical proportions; when measurement is impossible due to data limitations, alternative methods can be selected. k Let ∈[0.45,0.65] and let h = kL An approximate estimate is then made. Subsequently, the centroid penetration depth is obtained from the finite element simulation output. H p Then, press H i = H p - h Calculate the initial burial depth of the top of the mine H i and the H i As the starting burial depth for recursive superposition H 0 This ensures that the physical quantities superimposed on the subsequent sedimentation thickness sequence are consistent (both are "the burial depth of the mine top relative to the seabed surface").

[0088] Other technical solutions also include: S6. Conduct on-site reconnaissance and locate the number of [number missing] within the predicted sea area. N A suspected sea mine target. N For the stated integers, the integers are greater than 10; N Among the suspected sea mine targets K Excavation will proceed for each target. K greater than 1 and less than or equal to N The integer; the measurement obtained K The actual burial depth data of each target includes the measured burial depth value of each target; K The measured burial depth of each target is compared and analyzed with the predicted burial depth at the corresponding coordinate position in the final burial depth prediction distribution map.

[0089] To verify the accuracy of the prediction results in the above technical solution, the following field verification steps can be performed. After the prediction map is generated, a special survey operation is organized in the predicted sea area. Equipment such as multibeam echo sounders, side-scan sonar, magnetometers, or shallow seismic profilers are used to identify and locate suspicious seabed targets, forming a list of suspected targets and recording their geographic coordinates. Subsequently, based on the project schedule and risk assessment, targets are selected from the list. KVerification excavation was conducted on representative or key suspected targets. Overlying sediments were removed mechanically or hydraulically to expose the top of the target. The vertical distance from the current seabed surface to the top of the target was measured using a sounding rod or a high-precision distance sensor as the measured depth value. Finally, the measured points were spatially correlated with the predicted depth raster in GIS software. The predicted depth values ​​for the corresponding locations were extracted and compared pairwise with the measured values. Error statistics (such as mean error and root mean square error) were calculated to quantitatively assess the accuracy and reliability of the prediction model and provide a basis for subsequent model parameter correction and optimization.

[0090] Example 1 This embodiment takes the prediction of the burial depth of historical mines in the sea area of ​​a certain offshore wind power project on the east side of a strait as an example.

[0091] Step 1: Establish a three-dimensional finite element model of the seabed soil penetrated by the mine. 1.1 Determine the mine type and modeling parameters This embodiment targets the MK13 and MK26 bottom mines, which may have been historically deployed. Their numerical models are as follows: Figure 1 As shown.

[0092] 1.2 Creating the Geometric Model and Mesh A comprehensive 3D model containing mines and seabed soil, created in ABAQUS, is shown below. Figure 2 As shown in the figure, the seabed soil in this model is sufficiently sized to avoid boundary reflection interference.

[0093] Model mesh generation as follows Figure 3 As shown, the mines are represented by C3D10M elements, and the seabed soil is represented by EC3D8R Euler elements. The elements are further densified in the penetration influence zone, with a total of approximately 580,000 elements.

[0094] Step 2: Setting the material constitutive model and parameters 2.1 Selection of Constitutive Model for Seabed Soil The seabed surface soil in the study area is mainly composed of "clayous silt", and the Mohr-Coulomb elastoplastic constitutive model was adopted.

[0095] 2.2 Input Model Parameters Key parameters are shown in Table 1: Table 1. Parameters of the Finite Element Model of Seabed Soil Material <![CDATA[Density (kg / m 3 )]]> Angle of internal friction (°) Cohesion (kPa) clayey silt 1900 34 14 Step 3: Calculate the bottoming velocity and simulate the initial penetration depth 3.1 Calculate the bottom contact velocity at different water depths Based on the equations of motion of a mine entering the water, the contact velocities at different water depths are calculated as shown in Table 2. Figure 4 As shown.

[0096] Table 2 Vertical contact velocity of mines at different water depths Water depth (m) 0 5 10 20 30 40 Bottom contact velocity (m / s) 40 33 28 20 15 12 3.2 Perform CEL finite element analysis and obtain the initial penetration depth. Different bottoming velocities were used as initial conditions to input into the model for explicit dynamic analysis. Figure 5 The comparison of penetration depth and velocity attenuation of MK13 and MK26 mines under the same conditions shows that the MK13 mine penetrates deeper (1.64m vs 1.24m). Figure 6 The calculation model settings for different penetration angles (vertical and inclined) are shown.

[0097] The simulation results of the penetration depth are shown in Table 3 and Figure 7 As shown, the penetration depth decreases with increasing water depth, ranging from 0.64 m to 1.64 m.

[0098] Table 3. Penetration depth of mines at different water depths Water depth (m) 0 5 10 20 30 40 Penetration depth (m) 1.64 1.41 1.28 1.00 0.78 0.64 Step 4: Analyze the historical evolution of seabed erosion and deposition. 4.1 Acquiring and Processing Multi-Period Historical Nautical Charts Historical nautical charts from 1946, 1969, 1986, and 2012 were collected. After unifying the water depth to the theoretical depth datum, the underwater topographic distribution for each period was obtained as follows: Figures 8-11 .

[0099] 4.2 Calculation of scouring and silting variations and maximum siltation thickness Plotting changes in scour and sedimentation over different time periods ( Figures 12-14 ), and calculated the distribution of maximum sediment thickness between 1946 and 2012 ( Figure 15 The results showed that the siltation thickness was ≤4m in most areas, and up to 8m in some local areas.

[0100] 4.3 Identification and Handling of Disturbed Areas To quantify environmental disturbance factors, this embodiment identifies disturbed areas in the study area based on marine benthic organism survey data and seabed engineering records.

[0101] First, historical benthic organism survey data for this sea area were obtained. The data showed that the main disturbing species were Manila clams and smooth blue clams. Based on the biomass in the survey data, this embodiment quantifies and generates the marine biological activity intensity index B in the following specific manner: Parameter determination: Based on the ecological habits of the species, the disturbance weight W1=1.5 for *Mackerelia pennifolia* (burrowing depth >10cm) and the weight W2=1.0 for *Clam sclerotium* (burrowing depth 5-10cm). Survey data shows that the biomass of the two organisms in a representative area is A1=50 g / m³. 2 A2 = 30 g / m 2 .

[0102] Index Calculation: To normalize the index value, a standardization coefficient α = 0.01 is set. Calculated using the formula B = α·Σ(Wi·Ai), the value for this region is B = 0.01 × (1.5 × 50 + 1.0 × 30) = 0.01 × (75 + 30) = 1.05.

[0103] Spatialization and Threshold Setting: Similar calculations were performed for all survey sites, and a continuous "Bioactivity Intensity Index B Distribution Map" was generated using Kriging interpolation in ArcGIS software. A preset threshold B was set. t By statistically analyzing historical B values, this embodiment calculates B by taking the average B index of all sites plus 30%. t =0.80.

[0104] Secondly, historical records of seabed engineering activities were collected, and the scope (vector polygon) of the dredging projects carried out in the northwestern part of the study area between 1985 and 1995 was obtained.

[0105] Subsequently, the "Maximum Sedimentation Thickness Distribution Map" was analyzed in ArcGIS software. Figure 15 The disturbed area is identified, and the specific steps are as follows: Using the Raster Calculator, enter the conditional statement Con(“B-index distribution map”>0.80, 1, 0) to generate a binary raster of “Bio-disturbance zone”.

[0106] Convert the vector polygon of "Dredging Project Scope" into raster data and assign it the value 1.

[0107] Perform a Boolean OR operation on the two grids to obtain the final "Comprehensive Disturbance Area" label layer. Any grid that satisfies B>0.80 or is located within the historical project area is labeled as a "disturbance area".

[0108] Post-processing: In the final generated drawing (such as...) Figure 15 In the diagram, the marked area is covered with a light gray diagonal line that is different from the main image, and the legend is explained as follows: "The diagonal area is disturbed by biological or engineering activities, and the scouring and silting patterns may be abnormal."

[0109] Step 5: Coupled prediction of final mine burial depth First, based on the initial penetration depth H obtained from the simulation in step three...p (e.g., 1.64m) and the external dimensions of the mine determine the initial burial depth H. i For the MK13 type mine (total length) L =2.1m), axial distance from the center of mass to the top h Preferably, the value is determined by measurement using a three-dimensional geometric model; however, in this embodiment, due to the lack of a complete three-dimensional model suitable for measurement, an approximate value is used. h = kL Where k = 0.65, then h ≈ 1.37m, and thus H i = H p - h =1.64-1.37=0.27m.

[0110] Then, based on the water depth data from four periods (1946, 1969, 1986, and 2012), the thickness sequence of scour and deposition changes {Δ} was calculated. H 1,Δ H 2,Δ H 3}(correspond Figures 12-14 ),by H 0 = H i As the initial burial depth, according to H t = H t-1 +Δ H t By recursively calculating, the final burial depth distribution for the final year (2012) is obtained. H T At the same time, for { H 0, H 1, H 2, H 3 The maximum possible burial depth distribution within a historical period is obtained by taking the maximum value for each pixel. H max The final predicted distribution of mine burial depth is as follows: Figure 16 As shown in the figure. Predictions indicate that the burial depth in most areas is ≤3m, with some areas >4m and a maximum of approximately 7m.

[0111] 5.2 Dynamic Evolution Simulation If calculated recursively according to the sub-period scouring and silting sequence, the burial depth change process of the mine every ten years after its deployment can be simulated. Specifically, taking the initial burial depth H... i Starting from (0.27m), the values ​​are accumulated sequentially. Figure 12-14 The thickness of siltation changes at different time periods shown can be used to dynamically determine the possible burial depth of mines at any point in time, such as 1969 and 1986.

[0112] Step 6: Field Verification and Comparative Analysis Based on the predictions, scanning and excavation verification were carried out. A total of 60 suspected targets were found, of which 18 were excavated. Statistical analysis of the measured burial depth showed that 67% of the targets were <3m, only 2 were >4m, with a maximum depth of 6.8m. The measured distribution and... Figure 16 The high degree of consistency in the predicted trends verifies the reliability of this method.

[0113] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details and embodiments shown and described herein.

Claims

1. A method for predicting the burial depth of sea mines based on the CEL method and multi-year sedimentary evolution data, characterized in that, Includes the following steps: S1. Establish a three-dimensional finite element model of a mine penetrating the seabed soil. The model is constructed using the coupled Eulerian-Lagrange method, wherein the mine entity is simulated using a Lagrange mesh and the seabed soil is simulated using an Eulerian mesh. S2. Set the material constitutive model and corresponding parameters of the seabed soil in the three-dimensional finite element model, and set the contact and friction parameters between the mine and the soil; the parameters include the mine mass. M Frontal projection area S drag coefficient C y Drainage volume V dis Seawater density ρ w Soil density, internal friction angle, cohesion, and friction coefficient between the mine and the soil; S3. Based on the mine's mass, dimensions, drag coefficient, seawater density, discharge volume, and water depth, combined with the mine's initial entry velocity... v 0 and water entry angle θ Calculate the vertical component of the mine's bottom contact velocity for different water depths. V(X) The vertical component of the bottoming velocity was used as an initial condition and input into the three-dimensional finite element model for explicit dynamic analysis to simulate the penetration depth of the mine's center of mass in the seabed soil. H p ; S4. Obtain historical nautical chart depth data for the predicted sea area over a historical period, unify the historical nautical chart depth data to the same depth datum, and calculate the seabed erosion and deposition thickness sequence {Δ} for the predicted sea area over the historical period by comparing two adjacent periods. H t }; S5, Based on the centroid penetration depth H p The initial burial depth of the mine's top relative to the initial seabed surface is determined by the mine's external dimensions. H i ;by H 0 = H i As the initial burial depth, according to the aforementioned scouring and silting change thickness sequence {Δ H t By recursively superimposing the results, the final burial depth distribution for the termination year can be obtained. H T The maximum possible burial depth distribution within a historical time period was obtained by extracting the maximum value per pixel. H max Output a distribution map of the predicted burial depth.

2. The method as described in claim 1, characterized in that, The process of establishing the three-dimensional finite element model in step S1 includes: Based on the actual external dimensions of the target mine type, including its total length L With the maximum diameter D A three-dimensional geometric model of the mine entity was established and imported into the finite element analysis software; Based on mine quality M Assign material properties to the mine entity, including its material density. ρ mine According to the formula ρ mine = M / V model The calculation yielded, where V model Let V be the volume of the three-dimensional geometric model; An Eulerian grid domain of seabed soil is created around the mine entity, wherein the dimensions of the Eulerian grid domain in the longitudinal and transverse directions are not less than the maximum diameter of the mine. D 80 times the length of the mine, and its vertical dimension is not less than the total length of the mine. L 40 times; Define a general contact pair between a mine entity and seabed soil, and assign the friction coefficient between the mine and the soil to the tangential behavior property of the contact pair; Boundary conditions are applied to the Eulerian grid domain of the seabed soil. These boundary conditions include constraints on the normal displacements of all sides of the Eulerian grid domain and constraints on the bottom surface of the Eulerian grid domain. X Directional displacement, Y Directional displacement and Z Directional displacement.

3. The method as described in claim 1, characterized in that, In step S2, the material constitutive model of the seabed soil adopts the Mohr-Coulomb elastoplastic constitutive model; The parameters of the Mohr-Coulomb elastoplastic constitutive model include the density, internal friction angle, and cohesion of the seabed soil, with the density ranging from 1500 kg / m³. 3 Up to 2100 kg / m 3 The internal friction angle ranges from 25° to 40°, and the cohesion ranges from 5 kPa to 30 kPa.

4. The method as described in claim 1, characterized in that, In step S3, the vertical component of the mine's bottom-touching velocity for different water depths is calculated. V(X) The steps include: Based on the quality of the sea mine M Frontal projection area S drag coefficient C y Seawater density ρ w Drainage volume V dis and gravitational acceleration g Through formula Calculate the maximum speed of a naval mine in water. v np ; Based on the one-dimensional motion equation of the mine along its trajectory after entering the water, and combined with the initial entry velocity... v 0. Water entry angle θ With vertical water depth X Calculate the depth of the mine X velocity along the trajectory u(X) ;in, ； 'a' is the net acceleration constant caused by the difference between gravity and buoyancy when drag is ignored, and ; The velocity along the trajectory u(X) Projecting it vertically, we obtain the vertical component of the bottom-touching velocity. V(X) = u(X) cos θ .

5. The method as described in claim 1, characterized in that, In step S4, while acquiring historical nautical chart depth data for the predicted sea area over a historical period, marine benthic organism survey data for the predicted sea area are also acquired. Based on the biomass in the survey data, a marine biological activity intensity index is quantified and generated. B ; ； Where α is the standardization coefficient; W i For the first i Perturbation weight per unit biomass of benthic organisms; A i For the first i Biomass of benthic organisms per unit area.

6. The method as described in claim 5, characterized in that, Step S4 also involves acquiring records of seabed engineering activities in the predicted sea area; when the biological activity intensity index... B Greater than the preset threshold B t If there are records of submarine engineering construction, the corresponding area in the scouring and silting thickness data is marked as the disturbed area.

7. The method as described in claim 1, characterized in that, The specific process of recursive superposition in step S5 is as follows: based on the initial burial depth H i As the initial burial depth H 0 According to the set of thicknesses of scour and siltation changes {Δ H 1,Δ H 2,…,Δ H T Perform recursion: H t = H t-1 +Δ H t , ( t =1,2,…, T The burial depth distribution at the end of each sub-period was obtained. H t ; and according to H max =max{ H 0 , H 1 ,…, H T The maximum possible burial depth distribution within a historical time period is obtained.

8. The method as described in claim 1, characterized in that, The initial burial depth mentioned in step S5 H i Determined by the following formula: H i = H p - h ； in, h The distance is the axial distance from the center of mass of the mine to the top of the mine. h Determined by measurement using a three-dimensional geometric model of the mine; when measurement is not possible, take... h = kL ,in k ∈[0.45,0.65].

9. The method as described in claim 1, characterized in that, Also includes: S6. Conduct on-site detection and locate N suspected mine targets in the predicted sea area, where N is an integer greater than 10; excavate K targets among the N suspected mine targets, where K is an integer greater than 1 and less than or equal to N; measure the actual burial depth of the K targets. The actual burial depth value is compared and analyzed with the predicted burial depth value at the corresponding coordinate position in the burial depth prediction distribution map to assess the prediction error.

Images