A method for optimizing parameters of deep-shallow coupled surface model based on data assimilation

Abstract

The present invention provides a method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation, comprising: generating an initial surface unstructured grid; creating a surface deep-shallow coupling model input data set; constructing a tectonic grid feature map of a terrestrial evolution; inputting the initial surface unstructured grid, a paleoelevation data set, a precipitation data set and a tectonic grid feature map of a terrestrial evolution into the surface deep-shallow coupling model, and comparing a three-dimensional surface simulation result with a real environmental state in a geological history period to perform state estimation; constructing a vertical structure input file according to the three-dimensional surface simulation result and observation data in a geological history period, and correcting and updating a parameter set using the vertical structure input file; if the distribution of the corrected and updated parameter set is too narrow, perturbation processing is performed on the parameter set according to the sensitivity of the model to the parameters; and assimilation of all observation data is completed by continuously cyclically correcting and perturbation processing parameters.

Description

Technical Field

[0001] The present invention relates to the technical field of deep-shallow coupling model parameter optimization, and in particular to a method for optimizing surface deep-shallow coupling model parameters based on data assimilation. Background Art

[0002] The finite volume method, as the basic and key technology in fluid mechanics, plays a vital role in numerical algorithms. This method is based on the conservation equations in integral form and constructs an algorithm framework by controlling the physical conservation properties of the volume, rather than relying on traditional differential equations. It has demonstrated its strong applicability and flexibility in various application fields of fluid mechanics.

[0003] However, although the finite volume method has made remarkable achievements in the field of fluid mechanics, it still faces certain challenges in data processing, especially in the integration of geological data. The lack of necessary data constraints may prevent the method from fully realizing its potential in practical applications.

[0004] The deep-shallow coupled surface model, as an important tool for understanding and predicting geological evolution, plays a key role in integrating and assimilating observational data into the model. However, parameter uncertainty in the model often leads to deviations in the prediction results and affects the accuracy of the model. Traditional data assimilation methods mainly focus on the estimation of model states, while ignoring the direct optimization of model parameters. Summary of the invention

[0005] In view of the problems existing in the prior art, the present invention provides a method for optimizing the parameters of a surface deep-shallow coupling model based on data assimilation, which optimizes the coupling model parameters using the finite volume method and observation data, thereby improving the prediction accuracy of the deep-shallow coupling process in the geological model.

[0006] To achieve the above object, the present invention provides a method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation, comprising the following steps:

[0007] Step 1: Generate an initial surface unstructured grid;

[0008] Step 2: Create the input dataset of the surface deep-shallow coupled model, including the paleo-elevation input dataset and the precipitation input dataset;

[0009] Step 3, constructing a structural grid feature map of the surface variation;

[0010] Step 4: Input the initial surface unstructured grid, paleoelevation dataset, precipitation dataset and tectonic grid characteristic map of the earth's evolution into the surface deep-shallow coupling model, run the surface deep-shallow coupling model to obtain the surface three-dimensional simulation results, and compare the surface three-dimensional simulation results with the real environmental state in the geological history period to perform state estimation;

[0011] Step 5: construct a vertical structure input file based on the surface three-dimensional simulation results and the observation data of the geological history period, and use the vertical structure input file to calibrate and update the parameter set to reflect the influence of the observation data on the model parameters;

[0012] Step 6: If the distribution of the corrected and updated parameter set is too narrow, the parameter set is perturbed according to the sensitivity of the model to the parameters;

[0013] Step 7: Repeat steps 3 to 6 to simulate the evolution of the surface over time, and assimilate all observation data through continuous cycles of parameter correction and perturbation processing.

[0014] Optionally, the step 1 specifically includes:

[0015] Set the grid node coordinates, grid cell node index and adjacent index of each node, time interval for model running calculation and initial grid density described by Lagrangian. Based on this, construct and save the initial surface unstructured grid in the form of a triangular irregular network. Load the generated initial surface unstructured grid file into the corresponding software or use library functions for visualization and drawing to verify the validity of the grid spatial distribution.

[0016] Optionally, the step 2 specifically includes:

[0017] Set an initial value and initial standard deviation for each parameter in the deep-shallow coupled surface model;

[0018] Map the newly created Euler grid coordinates of the same resolution and projection mode to the coordinate system of the paleoelevation dataset and precipitation dataset, so that the coordinate system and spatial resolution of the paleoelevation dataset and precipitation dataset are consistent;

[0019] The Kriging interpolation method is used to interpolate the paleoelevation dataset and precipitation dataset on the Euler grid.

[0020] Optionally, the step 3 specifically includes:

[0021] Using GPlates software to reconstruct the plate and generate initial three-dimensional velocity field data, wherein the three-dimensional velocity field data includes velocity components of each plate in the X, Y, and Z directions in the global spherical coordinate system;

[0022] Mapping the initial three-dimensional velocity field data to each node of the initial surface unstructured grid, so that each node of the initial surface unstructured grid can obtain velocity values ​​in the X, Y, and Z directions;

[0023] The X, Y, and Z velocity values ​​of each grid node are multiplied by the specified time interval to calculate the displacement of each grid node within this time period, and the initial displacement data for constructing the grid feature map is obtained;

[0024] Mapping the displacement data from the global spherical coordinate system to the newly created Euler grid coordinate system;

[0025] The final construction mesh feature map is generated and saved using the displacement data in the transformed Euler mesh coordinate system.

[0026] Optionally, step 4 specifically includes:

[0027] The initial surface unstructured grid, paleoelevation dataset, precipitation dataset and tectonic grid feature map of the earth's evolution are averaged to obtain all the input conditions required to run the deep-shallow coupled surface model;

[0028] Run the surface forward model to simulate the evolution from the past geological history to the present, and obtain the simulation results of three-dimensional landforms, river networks, sediment fluxes, and delta evolution, wherein the surface forward model is a component of the surface deep-shallow coupling model;

[0029] The three-dimensional surface simulation results are compared with the actual environmental conditions during the geological history period to evaluate the simulation accuracy and evolution laws of the model.

[0030] Optionally, the step 5 specifically includes:

[0031] Calculate the difference between the surface 3D simulation results and the observation data of the geological history period, and create a vertical structure input file based on the difference;

[0032] The vertical displacement in the vertical structure input file is applied to the surface deep-shallow coupling model in stages and step by step using the linear interpolation method to calibrate the model parameter set.

[0033] Optionally, step 6 specifically includes:

[0034] Export the horizontal displacement data from GPlates software and apply it to the last run of the surface deep-shallow coupled model to perform horizontal displacement operations;

[0035] The entire computational region is divided into several discrete control volumes, and each control volume is used as a computational unit for model simulation;

[0036] The ancient elevation data are uniformly subjected to data parameter disturbance processing, and the specific formula is as follows:

[0037] H=H0+alpha*deltaH0,

[0038] Wherein, H is the vertical value of the structural grid feature map, H0 is the two-dimensional distribution value of the pre-estimated structural grid feature map, and deltaH0=0.1*H0 is assumed, where deltaH0 represents the disturbance value of the tectonic uplift of the geological variation, which is 10% of the initial tectonic uplift. Every 5 million years, H=H0+alpha*deltaH0 is substituted, where alpha is the disturbance coefficient, alpha=-1.0, -0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8, 1.0. A parallel simulation experiment is carried out for each disturbance coefficient, and a total of 10 parallel numerical simulation experiments are carried out.

[0039] The vertical values ​​of the 10 5-million-year tectonic grid feature maps and the elevation values ​​of each point of the corresponding simulation results were statistically analyzed to obtain the Pearson correlation coefficient K between the two, and the intercept b of the regression line was calculated. The specific formula is as follows:

[0040] ,

[0041] ,

[0042] in, x is the vertical value of 10 5 million-year tectonic grid feature maps, y is the elevation value of each point in the simulation result, and They are x and y The sample mean of ;

[0043] Calculate the difference between the simulation results calculated when the perturbation coefficient alpha is 0 and the simulation results calculated when the perturbation coefficient alpha is not 0 Y , and according to the formula Y =K* X + b Calculate △ X ;

[0044] △ X Applied to the vertical structure input file of the next time step, used for running the surface forward model simulation at the next time step.

[0045] After adopting the above technical solution, the present invention has at least the following beneficial effects:

[0046] 1. Improved the accuracy of parameter estimation of the deep-shallow coupled surface model

[0047] The present invention performs parameter perturbations by assuming true values ​​and collects a large number of observations, and performs data assimilation after comparing the differences between the two, and performs grid redivision processing at the same time. This improvement integrates paleotopographic maps, paleoprecipitation grids and horizontal displacement files into the same process, which significantly improves the accuracy of parameter estimation. Since the parameter perturbation and data assimilation processes effectively reduce the difference between the model and the real geological data, it ensures that the simulation results are closer to the actual geological evolution process. Therefore, a more accurate parameter estimation effect is produced, and the simulation ability of geological structure evolution is enhanced.

[0048] 2. Improved the accuracy and continuity of the deep-shallow coupled surface model simulation

[0049] The present invention develops a more reliable model set by optimizing the disturbance parameters, thereby improving the accuracy of parameter estimation in the simulation. This improvement not only improves the accuracy of data assimilation, but also effectively improves the continuity of surface sedimentary materials transported by rivers. Because of this coupling process that continuously corrects the sedimentation and erosion processes, the forward model is more consistent with the evolution process of the real geological history, providing more reliable and detailed surface process simulation results, enabling the model to more accurately reflect the details of the surface tectonic evolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0051] Figure 1 A schematic flow chart of a method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation provided in an embodiment of the present disclosure;

[0052] Figure 2 is the error diagram between the simulated paleoelevation and the observed paleoelevation after multiple data assimilation;

[0053] Figure 3 This is a three-dimensional visualization of the surface simulation results of Taiwan and its neighboring areas 3 million years ago. DETAILED DESCRIPTION

[0054] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0055] In the field of landform evolution simulation, the finite volume method has become the preferred numerical method. Before calculation, this method divides the calculation area into a series of non-repeating control volumes, each of which is defined by a representative node. By integrating the space and time of the conservation differential equation in the control volume and a specific time interval, the physical phenomenon can be accurately simulated. The unknown values ​​are at the grid points. By assuming the law of their change between the grid points, the function and its derivatives are interpolated in time and space or the assumption of the change line is made, and a set of discrete equations are obtained to integrate the control volume. The finite volume method constructs discrete equations from a physical point of view to ensure that each equation reflects the conservation of a certain physical quantity on the finite volume. The derivation process of this method is clear, the coefficients of the discrete equations have clear physical meanings, and ensure that the equations have conservation characteristics. The discrete equations of the finite volume method require the integral conservation of the dependent variable, which is applicable to any control volume and ensures the conservation of the entire calculation area. This is its significant advantage, and its more flexible assumptions can overcome the shortcomings of Taylor expansion discretization and have good adaptability to the grid. Especially when conducting surface fluid-solid coupling analysis, it can solve complex surface simulation problems well.

[0056] However, although the finite volume method has made remarkable achievements in the field of fluid mechanics, it still faces certain challenges in multi-layer coupling and data processing, especially in the integration of geological data. The lack of necessary data constraints may prevent the method from fully realizing its potential in practical applications.

[0057] Data assimilation is a mathematical method that combines observational data and numerical models to estimate the state of a system. This method has been widely used in meteorological forecasting, oceanography, environmental science, geophysics and other fields. The mathematical principles of data assimilation are mainly based on probability theory and mathematical statistics, involving mathematical concepts such as random variables, Bayesian inference, and optimization theory. The observational data and model predictions in data assimilation can be regarded as random variables. Random variables are quantities with uncertainty, and their probability distribution describes the uncertainty of possible values. In data assimilation, observational data usually contains observation errors, while model predictions contain model errors. These errors can be described by random variables. From an algorithmic perspective, data assimilation and artificial intelligence have a strong commonality, that is, how to make reasonable guesses. There are two main types of guesses: the first is to guess the relationship between model parameters and simulation results; the second is to estimate the specific values ​​of parameters.

[0058] The inventors have deeply explored the parameterization characteristics of the deep and shallow coupled surface model, and adopted a novel parameter optimization strategy for the specific sources of model deviations. By initializing the parameter set and performing state estimation, it is ensured that the uncertainty of the model state is fully constrained by the observation data. After the state reaches quasi-equilibrium, the difference between the parameter set and the observation set is calculated to update the parameter set. In the case where the distribution of the parameter set is too narrow, the parameter set is perturbated to maintain its adjustability. The accuracy of the parameters is optimized by continuously cycling parameter correction and perturbation processing until the assimilation of all observation data is completed.

[0059] The present invention is described in detail with specific embodiments below:

[0060] refer to Figure 1 The present disclosure provides a method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation, comprising the following steps:

[0061] Step 1: Generate an initial surface unstructured grid;

[0062] In this step, the key information to be prepared includes: grid node coordinates, which are used to define the specific position of each grid node in space; grid cell node index, which is used to identify the set of nodes that constitute each cell (Cell); each node adjacent index, which is used to clarify the relationship between each node and its adjacent nodes; setting the time interval for model running calculations, in units of million years (million years); setting the initial grid density described by Lagrangian, and based on this, constructing and saving the initial surface unstructured grid in the form of a triangular irregular network to ensure that the data better fits the surface morphology and is applied to the subsequent state estimation step. Before proceeding to the subsequent steps, the generated initial surface unstructured grid file is loaded into the corresponding software or the library function is used for visual drawing to verify the validity of the grid spatial distribution, ensure the accuracy and reliability of the grid data, and lay the foundation for subsequent simulation and analysis.

[0063] Step 2: Create the input dataset of the surface deep-shallow coupled model, including the paleo-elevation input dataset and the precipitation input dataset;

[0064] In this step, an initial value and initial standard deviation are first set for each parameter in the deep-shallow coupled surface model to ensure that the parameter set has an appropriate starting point and uncertainty range;

[0065] Then, when processing paleoelevation datasets (paleoelevation data from the geological history period, such as the altitude of a certain location 1 million years ago) and precipitation datasets, we first need to note that they have different coordinate systems. Paleoelevation datasets are generally grid data in the Mercator projection coordinate system, and paleoprecipitation data use the global spherical coordinate system. However, the data structures of these two datasets are similar. Paleoelevation datasets contain longitude, latitude, and paleoelevation; precipitation datasets also contain longitude, latitude, and precipitation. In view of this particularity, the newly created Euler grid coordinates of the same resolution and projection method are mapped to the coordinate system of the paleoelevation dataset and precipitation dataset, so that the coordinate system and spatial resolution of the paleoelevation dataset and precipitation dataset are consistent. This step is the key to identifying and utilizing the characteristics of different datasets.

[0066] The Kriging interpolation method is used to interpolate the paleo-elevation dataset and precipitation dataset on the Euler grid to ensure that the information extracted from the original dataset can be accurately represented on the new grid. This step is a key step in extracting and utilizing the characteristics of the paleo-dataset and provides important data support for subsequent geological simulation and analysis.

[0067] Step 3, constructing a structural grid feature map of the surface variation;

[0068] In this step, GPlates software is first used to reconstruct the plate and generate initial three-dimensional velocity field data. The three-dimensional velocity field data includes the velocity components of each plate in the X, Y, and Z directions in the global spherical coordinate system (unit: cm / year). These data are provided in the form of discrete point data, and each point corresponds to the velocity information of a geographical location. The velocity field data will provide the necessary original velocity information for the next step of allocating velocity values ​​on the structural grid;

[0069] The initial three-dimensional velocity field data is mapped to each node of the initial surface unstructured grid. To this end, a spatial function is introduced to establish the spatial correspondence between each grid node and the GPlates three-dimensional velocity field data point. Specifically, the nearest neighbor interpolation method is used to find the nearest neighbor node on the construction grid for each velocity data point provided by GPlates according to its geographical location, and the velocity value is assigned to this node, that is, each velocity data point provides a direct velocity source for the construction grid, so that each node of the initial surface unstructured grid can obtain the velocity value in the X, Y, and Z directions;

[0070] After obtaining the X, Y, and Z velocity data of each grid node, the X, Y, and Z velocity values ​​of each grid node are multiplied by the specified time interval to calculate the displacement of each grid node within this time period, which is equivalent to converting the instantaneous velocity of each node into a long-term displacement, that is, the tectonic uplift of each node, and obtaining the initial displacement data of the tectonic grid feature map. The velocity data in the previous step is used here to determine the movement of each node on the geological time scale;

[0071] The displacement data generated in the previous step are still in the global spherical coordinate system (i.e., expressed in latitude, longitude and spherical coordinates). In order to better express these data in the construction grid feature map, the displacement data need to be mapped from the global spherical coordinate system to the newly created Euler grid coordinate system. Through coordinate conversion, the displacement data in these spherical coordinate systems are recalculated into the Euler grid coordinate system, so that the entire construction process can be expressed in the standard framework of the Euler grid coordinate system. This step ensures that the surface motion data can be presented in a coordinate system that is more suitable for scientific analysis, and ensures that the displacement data remains accurate in the new coordinate system.

[0072] The final tectonic grid feature map is generated and saved using the displacement data in the converted Euler grid coordinate system. The tectonic grid feature map summarizes the displacement information of all nodes and fully describes the evolution of surface tectonic movement. In order to check whether the data processing is correct, the displacement velocity file should be output as a visualization file and these tectonic grid feature maps should be stored.

[0073] Through the above steps, the generation of the three-dimensional velocity field and the final creation of the grid feature map are completed.

[0074] Step 4: Input the initial surface unstructured grid, paleoelevation dataset, precipitation dataset and tectonic grid characteristic map of the earth's evolution into the surface deep-shallow coupling model, run the surface deep-shallow coupling model to obtain the surface three-dimensional simulation results, and compare the surface three-dimensional simulation results with the real environmental state in the geological history period to perform state estimation;

[0075] In the million-year time scale of surface tectonic simulation, one of the core ideas is to average the model input data, that is, to calculate the annual mean. In the actual geological evolution process, the erosion rate, precipitation, etc. fluctuate in different months and seasons of the year, so the changes in 12 months are averaged to obtain the average state for the whole year. After this processing, although it cannot fully reflect every moment of change in the geological system, the average value can more effectively smooth these changes and capture the long-term trend of the geological evolution process. This method has advantages in ensuring the stability and operability of the calculation, and at the same time provides reasonable initial input data for long-term simulation.

[0076] In this step, the initial surface unstructured grid, paleoelevation dataset, precipitation dataset, and tectonic grid feature map of the earth's evolution are first averaged to obtain all the input conditions required to run the deep-shallow coupled surface model;

[0077] Run the surface forward model to simulate the evolution from the past geological history to the present, and obtain the simulation results of three-dimensional landforms, river networks, sediment fluxes, and delta evolution, wherein the surface forward model is a component of the surface deep-shallow coupling model;

[0078] At the same time, the surface forward model does not operate independently. It is combined with the data assimilation process to form a coupled model. During the operation of the coupled model, the surface forward model continuously receives feedback from data assimilation, compares the three-dimensional surface simulation results with the actual environmental conditions during the geological history period, and evaluates the simulation accuracy and evolution laws of the model. This process not only ensures the accuracy of the forward model, but also provides a more reliable prediction capability for the surface evolution process.

[0079] Step 5: construct a vertical structure input file based on the surface three-dimensional simulation results and the observation data of the geological history period, and use the vertical structure input file to calibrate and update the parameter set to reflect the influence of the observation data on the model parameters;

[0080] Over the past few decades, with the continuous deepening of paleoelevation research, a large amount of paleoelevation and paleo-latitude and longitude data about the geological history period has been accumulated, which provides an important calibration basis for the simulation results. In order to evaluate the differences between the model simulation results and the real paleoelevation data in different periods, the simulated three-dimensional surface simulation results are first obtained and compared with the observation data of the geological history period. Through the difference between the two, a vertical displacement paleoelevation map is created (the corresponding model data file is also called a vertical structure input file). It is assumed that this vertical displacement paleoelevation map has taken into account the steady-state response of dynamic terrain changes, crustal bending and other tectonic actions, reflecting the influence of geological structure on the evolution of the earth's surface. In the geological history period, the sedimentation and erosion processes directly affect the accuracy of paleoelevation simulation. Therefore, the paleosurface sedimentary materials in the three-dimensional simulation will be transported by rivers, and the finite volume method is used for calculation to ensure the accurate simulation of the sedimentation and erosion processes. However, in order to more accurately reflect the impact of vertical structure, we will not directly apply all vertical displacements in the vertical structure input file to the simulation at one time, because these displacement values ​​represent evolutionary changes over a long time scale (millions of years). If they are applied at one time, although the computational efficiency will be improved, it may lead to obvious deviations and inaccuracies in simulation results such as river sediment flux.

[0081] To this end, a linear interpolation method is used to apply the vertical displacement in stages and step by step. For example, for simulations between 20 and 19 million years, a vertical displacement of 0.2 is applied; for simulations between 20 and 18 million years, a vertical displacement of 0.4 is applied, and so on. In this way, the simulation can gradually and smoothly reflect the historical evolution of vertical structures in each time period. In addition, this vertical displacement coefficient can be flexibly adjusted manually in time and space according to the nonlinear characteristics of the tectonic history to ensure that the model more accurately reflects the real change process of geological structures.

[0082] Step 6: If the distribution of the corrected and updated parameter set is too narrow, the parameter set is perturbed according to the sensitivity of the model to the parameters;

[0083] Export the horizontal displacement data from the GPlates software and apply it to the last run of the deep-shallow coupled surface model. Save the mesh topology, geometry, and attribute information of the surface model to a file. This information includes the element type, node coordinates, longitude and latitude of the mesh, and horizontal displacement data. After performing the horizontal displacement operation, the surface model needs to be re-meshed to ensure that it matches the new displacement data.

[0084] Although we are dealing with a surface model, surface process simulation is actually a three-dimensional computational process, because the spatial parameters of the surface are not only composed of longitude and latitude coordinates, but also paleo-elevation values. Therefore, when re-meshing, the entire computational area needs to be divided into discrete control volumes, usually in the form of a structured grid. The computational area refers to the entire geological area involved in the surface process simulation, which is divided into many small control volumes for processing.

[0085] The control volume is the basic unit of numerical calculation. It represents a closed space area and can be presented as a cube, hexahedron, tetrahedron or other geometric shapes. In this way, the control volume, as the smallest calculation unit of the simulation, can more accurately capture the spatial variation characteristics of the surface. Generally, the finer the division, the higher the accuracy of the model simulation. Re-dividing the grid and using discrete control volumes can more effectively fit and reproduce the paleo-geomorphic characteristics, thereby improving the accuracy of the simulation results.

[0086] In order to deeply explore the parameterization characteristics of the deep-shallow coupled model of the surface, a novel parameter optimization strategy was adopted for the specific source of model deviation, that is, paleoelevation data, and the paleoelevation data was uniformly perturbed. The specific formula is as follows:

[0087] H=H0+alpha*deltaH0,

[0088] Wherein, H is the vertical value of the structural grid feature map, H0 is the two-dimensional distribution value of the pre-estimated structural grid feature map, and deltaH0=0.1*H0 is assumed, where deltaH0 represents the disturbance value of the tectonic uplift of the geological variation, which is 10% of the initial tectonic uplift. Every 5 million years, H=H0+alpha*deltaH0 is substituted, where alpha is the disturbance coefficient, alpha=-1.0, -0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8, 1.0. A parallel simulation experiment is carried out for each disturbance coefficient, and a total of 10 parallel numerical simulation experiments are carried out.

[0089] For the vertical values ​​of the 10 5-million-year tectonic grid characteristic maps, it is a disturbance state, which is the independent variable. The elevation value of each point of the corresponding simulation result, that is, the corresponding dependent variable, is obtained through statistical analysis. The Pearson correlation coefficient K of the two is obtained, and the intercept b of the regression line is calculated. The specific formula is as follows:

[0090] ,

[0091] ,

[0092] in, x is the vertical value of 10 5 million-year tectonic grid feature maps, y is the elevation value of each point in the simulation result, that is, the corresponding dependent variable, and They are x and y The sample mean of ;

[0093] Calculate the difference between the simulation results calculated when the perturbation coefficient alpha is 0 and the simulation results calculated when the perturbation coefficient alpha is not 0 Y , and according to the formula Y =K* X + b Calculate △ X ;

[0094] In the surface simulation process, in order to improve the consistency between the simulation results and the actual geological history, we distribute the difference between the simulated elevation calculated at each 10,000-year interval and the real elevation calibration point to each grid point of the tectonic grid, and reapply this calibrated difference (△X calculated in the previous step) to the vertical structural input file of the next time step, which is directly used for the surface forward model operation simulation of the next time step. In this way, the continuity of the ancient surface sedimentary materials transported by rivers in the surface simulation has been significantly improved, so that the simulation better reflects the real process of river transportation and deposition of materials.

[0095] Step 7: Repeat steps 3 to 6 to simulate the evolution of the surface over time, and assimilate all observation data through continuous cycles of parameter correction and perturbation processing.

[0096] Repeat steps 3 to 6, and continuously assimilate data through this iterative process to simulate the evolution of surface structures over time, and complete the assimilation of all observation data through continuous cycles of parameter correction and disturbance processing. This not only improves the simulation accuracy of geological history processes, but also optimizes the simulation results.

[0097] The method described in the above embodiment is used to simulate the geomorphic changes in Taiwan and its neighboring areas 3 million years ago, and the simulation results of the geomorphic changes in Taiwan and its neighboring areas 3 million years ago are obtained. Figure 2 The error between the simulated paleo-elevation after multiple data assimilation and the observed paleo-elevation is shown. Figure 3 The results of the three-dimensional visualization of the surface simulation of Taiwan and its neighboring areas 3 million years ago are shown. Figure 2 It can be seen that the error range between the simulated paleoelevation and the observed paleoelevation is mainly concentrated between -10 meters and 10 meters, which shows that the simulation accuracy of the model has been significantly improved, which can effectively help scientific researchers to conduct refined marine geology and paleocean research. Based on this data assimilation model, it is observed that the uplift of the Taiwan orogenic belt at 5Ma strengthened the transport of sediments by the (paleo) Taiwan Warm Current, and promoted the settlement of sediments in the Taixi Basin and the Southwest Taixi Basin, providing important clues for understanding the geological evolution of the region.

[0098] Although the present invention has been disclosed as above by way of embodiments, it is not intended to limit the present invention. Any person skilled in the art may make some changes and modifications without departing from the spirit and scope of the present invention. Therefore, the protection scope of the present invention shall be determined by the claims.

Claims

1. A method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation, characterized in that: The following steps are involved: Step 1: Generate an initial surface unstructured grid; Step 2: Create the input dataset of the surface deep-shallow coupled model, including the paleo-elevation input dataset and the precipitation input dataset; Step 3, constructing a structural grid feature map of the surface variation; Step 4: Input the initial surface unstructured grid, paleoelevation dataset, precipitation dataset and tectonic grid characteristic map of the earth's evolution into the surface deep-shallow coupling model, run the surface deep-shallow coupling model to obtain the surface three-dimensional simulation results, and compare the surface three-dimensional simulation results with the real environmental state in the geological history period to perform state estimation; Step 5: construct a vertical structure input file based on the surface three-dimensional simulation results and the observation data of the geological history period, and use the vertical structure input file to calibrate and update the parameter set to reflect the influence of the observation data on the model parameters; Step 6: If the distribution of the corrected and updated parameter set is too narrow, the parameter set is perturbed according to the sensitivity of the model to the parameters; Step 7: Repeat steps 3 to 6 to simulate the evolution of the surface over time, and assimilate all observation data by continuously cyclically correcting parameters and performing disturbance processing; The step 6 specifically includes: Export the horizontal displacement data from GPlates software and apply it to the last run of the surface deep-shallow coupled model to perform horizontal displacement operations; The entire computational region is divided into several discrete control volumes, and each control volume is used as a computational unit for model simulation; The ancient elevation data are uniformly subjected to data parameter disturbance processing, and the specific formula is as follows: H=H0+alpha*deltaH0, Wherein, H is the vertical value of the structural grid feature map, H0 is the two-dimensional distribution value of the pre-estimated structural grid feature map, and deltaH0=0.1*H0 is assumed, where deltaH0 represents the disturbance value of the tectonic uplift of the geological variation, which is 10% of the initial tectonic uplift. H=H0+alpha*deltaH0 is substituted every 5 million years, where alpha is the disturbance coefficient, alpha=-1.0, -0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8, 1.

0. A parallel simulation experiment is carried out for each disturbance coefficient, and a total of 10 parallel numerical simulation experiments are carried out. The vertical values ​​of the 10 5-million-year tectonic grid feature maps and the elevation values ​​of each point of the corresponding simulation results were statistically analyzed to obtain the Pearson correlation coefficient K between the two, and the intercept b of the regression line was calculated. The specific formula is as follows: Where x is the vertical value of the tectonic grid characteristic map of 10 times 5 million years, and y is the elevation value of each point in the simulation result. and are the sample means of x and y respectively; Calculate the difference △Y between the simulation result calculated when the disturbance coefficient alpha is 0 and the simulation result calculated when the disturbance coefficient alpha is not 0, and calculate △X according to the formula Y=K*X+b; Apply △X to the vertical structure input file of the next time step and use it to run the simulation of the surface forward model of the next time step.

2. The method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation according to claim 1, characterized in that: The step 1 specifically includes: Set the grid node coordinates, grid cell node index and adjacent index of each node, time interval for model running calculation and initial grid density described by Lagrangian. Based on this, construct and save the initial surface unstructured grid in the form of a triangular irregular network. Load the generated initial surface unstructured grid file into the corresponding software or use library functions for visualization and drawing to verify the validity of the grid spatial distribution.

3. The method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation according to claim 1 is characterized in that: The step 2 specifically includes: Set an initial value and initial standard deviation for each parameter in the deep-shallow coupled surface model; Map the newly created Euler grid coordinates of the same resolution and projection mode to the coordinate system of the paleoelevation dataset and precipitation dataset, so that the coordinate system and spatial resolution of the paleoelevation dataset and precipitation dataset are consistent; The Kriging interpolation method is used to interpolate the paleoelevation dataset and precipitation dataset on the Euler grid.

4. The method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation according to claim 1, characterized in that: The step 3 specifically includes: Using GPlates software to reconstruct the plate and generate initial three-dimensional velocity field data, wherein the three-dimensional velocity field data includes velocity components of each plate in the X, Y, and Z directions in the global spherical coordinate system; Mapping the initial three-dimensional velocity field data to each node of the initial surface unstructured grid, so that each node of the initial surface unstructured grid can obtain velocity values ​​in the X, Y, and Z directions; The X, Y, and Z velocity values ​​of each grid node are multiplied by the specified time interval to calculate the displacement of each grid node within this time period, and the initial displacement data for constructing the grid feature map is obtained; Mapping the displacement data from the global spherical coordinate system to the newly created Euler grid coordinate system; The final construction mesh feature map is generated and saved using the displacement data in the transformed Euler mesh coordinate system.

5. The method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation according to claim 1, characterized in that: The step 4 specifically includes: The initial surface unstructured grid, paleoelevation dataset, precipitation dataset and tectonic grid feature map of the earth's evolution are averaged to obtain all the input conditions required to run the deep-shallow coupled surface model; Run the surface forward model to simulate the evolution from the past geological history to the present, and obtain the simulation results of three-dimensional landforms, river networks, sediment fluxes, and delta evolution, wherein the surface forward model is a component of the surface deep-shallow coupling model; The three-dimensional surface simulation results are compared with the actual environmental conditions during the geological history period to evaluate the simulation accuracy and evolution laws of the model.

6. The method for optimizing parameters of a surface deep-shallow coupling model based on data assimilation according to claim 1, characterized in that: The step 5 specifically includes: Calculate the difference between the surface 3D simulation results and the observation data of the geological history period, and create a vertical structure input file based on the difference; The vertical displacement in the vertical structure input file is applied to the surface deep-shallow coupling model in stages and step by step using the linear interpolation method to calibrate the model parameter set.

Images