A quadrilateral mesh optimization method for geological models based on particle swarm optimization algorithm

Abstract

The invention discloses a method for optimizing a quadrilateral mesh of a geological model based on a particle swarm optimization algorithm, and relates to the technical field of earthquake numerical simulation, including: reading the quadrilateral mesh data of the geological model, and storing it as a mesh node table and a quadrilateral table; adding the nodes to be optimized to the node table to be optimized; determining the solution space for all the mesh nodes in the node table to be optimized, and initializing the initial position, individual optimal position, and global optimal position of all particles in the particle swarm in the solution space; randomly updating the particle velocity for all particles in the particle swarm, and balancing the degree of approach of the particles to the individual optimal position and the global optimal position based on the weight parameter; after a finite number of optimization iterations, determining the global optimal position as the position after the node optimization; and saving the optimized quadrilateral mesh data of the geological model. The method can solve the problems of low calculation accuracy, low efficiency, and poor stability caused by the quadrilaterals with poor shapes in the finite element calculation of the geological model.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake numerical simulation, and in particular to a method for optimizing quadrilateral grids of geological models based on a particle swarm optimization algorithm. Background Art

[0002] At present, with the increase of engineering needs and application scenarios, structural finite element calculation has been more widely used in dealing with earthquake numerical simulation. In the process of earthquake research, it is necessary to obtain the cross-section of the seismic zone, and then mesh the cross-section to obtain a quadrilateral mesh model. The quadrilateral mesh model can help simulate the propagation of seismic waves under these complex conditions, thereby predicting possible earthquake response characteristics and providing guidance for subsequent data processing and interpretation. At present, when meshing the cross-section of the seismic zone, the quadrilateral mesh generated by the traditional meshing algorithm often has poorly shaped quadrilaterals when processing complex models. These quadrilaterals have a certain impact on the accuracy, efficiency and stability of the calculation. Optimizing the deformed quadrilaterals by moving the quadrilateral mesh nodes can obtain a better finite element model and provide high-quality data support for subsequent finite element calculations. The existing means of optimizing quadrilateral meshes are mostly to optimize the mesh generation algorithm, adjust the mesh parameters or adopt a better meshing technology. The degree of optimization of the quadrilateral mesh is very limited and cannot fully meet the actual needs. Therefore, it is urgent to design an efficient, low-cost and intelligent method for optimizing the shape of quadrilateral grids in geological models. Summary of the invention

[0003] The present invention provides a method for optimizing quadrilateral meshes of geological models based on a particle swarm optimization algorithm, which can solve the problems of low calculation accuracy, low efficiency and poor stability caused by poorly shaped quadrilaterals in finite element calculations of geological models.

[0004] In order to alleviate the above-mentioned problems, the technical solution adopted by the present invention is as follows:

[0005] The present invention provides a method for optimizing a quadrilateral mesh of a geological model based on a particle swarm optimization algorithm, comprising the following steps:

[0006] S101, read the geological model grid file obtained by the dissection, and store the read file data in an internal data structure composed of a grid node table and a quadrilateral table;

[0007] S102, traverse all nodes in the grid node table except the boundary nodes, add the nodes to be optimized to the node table to be optimized L, if the length of the node table to be optimized L is 0, jump to step S113, otherwise continue to execute step S103;

[0008] S103, start traversing the node table L to be optimized, and initialize the node number i=1;

[0009] S104, determine whether i is greater than the length of the node table L to be optimized, if so, jump to step S102 and perform the next round of optimization on the entire grid, otherwise continue to step S105;

[0010] S105, select the i-th node N from the node list L to be optimized i , determine the node N i The solution space of the particle swarm is initialized in the solution space, and the loss value of the initial position is calculated according to the loss function of the optimization target. Then, the individual optimal position and the global optimal position are initialized, and it is ensured that all particles will not leave the solution space during the movement.

[0011] S106, start iterative optimization of node N i Position, initialize the iteration round j = 1;

[0012] S107, determine whether the iteration round j is greater than the set maximum optimization iteration number M, if so, assign a value to i, i=i+1, and jump to step S104, otherwise continue to execute step S108;

[0013] S108, in the jth iteration, the particle number is initialized to k = 1, and the loss value of the optimal position of the individual particle is LocalLoss k ;

[0014] S109, determine whether the particle number k is greater than the set total number of particles Q, if so, assign a value to j, j=j+1, and jump to step S107, otherwise continue to execute step S110;

[0015] S110, randomly generate particles P k The speed of the particle is calculated and the particle position is updated (to ensure that the particle position is confined to the solution space). During the position update process, the particle's approach to the individual optimal position and the global optimal position is balanced based on the weight parameter, and the particle P is calculated. k The loss value Loss at the latest position;

[0016] S111, if the loss value Loss at the latest position is less than the loss value LocalLoss at the optimal position of the particle individual k , then update the current position to the individual optimal position of particle k, and continue to execute step S112, otherwise the kth particle is optimized, assign a value to k, k=k+1, and return to step S109;

[0017] S112, if the loss value Loss at the latest position is less than the loss value at the global optimal position of the particle swarm, then update the current position to the global optimal position of the particle swarm, the kth particle is optimized, assign a value to k, k=k+1, and return to step S109;

[0018] S113, saving the optimized quadrilateral grid data of the geological model.

[0019] In a preferred embodiment of the present invention, step S101 specifically includes:

[0020] Parse and load mesh information from the standard mesh file used by the seismic finite element model. The mesh information includes the numbers and coordinates of mesh nodes, the numbers of quadrilaterals and their constituent nodes, and the numbers of mesh nodes on the boundaries of the seismic finite element model.

[0021] Traverse all grid node information, and store three types of information, namely, the number of the grid node, the coordinates of the grid node, and whether the grid node is at the boundary of the seismic finite element model, in the grid node table;

[0022] All quadrilateral information is traversed, and the quadrilateral number and the numbers of the four constituent nodes are stored in the quadrilateral table, and the quadrilateral number is stored in the information of the four constituent nodes in the grid node table.

[0023] In a preferred embodiment of the present invention, in step S102, the nodes to be optimized include four categories, namely a, b, c and d. The four categories of nodes are nodes that cause topological errors, nodes that cause concave quadrilaterals, nodes that cause quadrilaterals to degenerate into triangles, and nodes that cause angles that are too large or too small. The optimization order priority of the four categories of nodes to be optimized decreases from a to d.

[0024] In a preferred embodiment of the present invention, the method for determining whether a certain node causes a topological error includes:

[0025] Find all quadrilaterals with the node as a vertex, and find all line segments on the quadrilateral that are not connected to the node. All line segments can enclose a closed area.

[0026] Take this node as the starting point and draw a ray in the negative direction of the y axis;

[0027] Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node has a topological error.

[0028] In a preferred embodiment of the present invention, the method for determining whether a node forms a concave quadrilateral includes:

[0029] Find all nodes connected to this node, and connect these nodes to form a closed area;

[0030] Take this node as the starting point and draw a ray in the negative direction of the y axis;

[0031] Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node causes the appearance of a concave quadrilateral.

[0032] In a preferred embodiment of the present invention, the node N is calculated i Ways to lose value at a position include:

[0033] Step 1: Calculate the number of nodes N i is the loss value of all quadrilaterals at the vertex, and the formula is as follows:

[0034]

[0035] Among them, loss quad is the loss value of a single quadrilateral; TopologyError indicates that there is a topological error in the quadrilateral, and two of the four sides are crossed. TopologyErrorLoss is the loss value of the quadrilateral in this case; ConcaveError indicates that the quadrilateral is a concave quadrilateral, and ConcaveErrorLoss is the loss value of the quadrilateral in this case; AngleError indicates that the vertex angle of the quadrilateral is too large or too small, and AngleErrorLoss is the loss value of the quadrilateral in this case; angle i is the i-th interior angle of the quadrilateral; len i is the ith ordered edge of the quadrilateral, It represents the ratio of two opposite sides in a quadrilateral;

[0036] Step 2: Calculate the number of nodes N i The loss value is as follows:

[0037]

[0038] Among them, loss node is the loss value of the node, n is the number of quadrilaterals with the node as the vertex, is the loss value of the ith quadrilateral, m is the number of edges with this node as the endpoint, len j and len k are the lengths of the j-th and k-th edges.

[0039] In a preferred embodiment of the present invention, in step S105, for node N i , find all nodes N i Connected nodes, these nodes are connected to form a closed area, that is, node N isolution space; calculate the loss value of the initial position of each particle according to the loss function of the optimization objective, and initialize the initial position with the lowest loss value as the global optimal position; for a certain particle, take the initial position of the particle as the initial individual optimal position of the particle.

[0040] In a preferred embodiment of the present invention, for a particle p k , the optimization process includes:

[0041] Calculate particle p k The loss value Loss at the latest position;

[0042] If the loss value Loss is greater than or equal to its individual optimal position, then particle p k Complete the optimization, otherwise update the current position to particle p k The individual optimal position of

[0043] If the loss value Loss is greater than or equal to the global optimal position, then particle p k Complete the optimization, otherwise update the current position to the global optimal position and complete the particle p k Optimization.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] (1) The finite element method uses interpolation functions to approximate the physical quantities within the unit. For regular-shaped quadrilateral units, interpolation functions can be constructed more easily and can more accurately approximate the true solution. By using a mesh generation method based on a particle swarm optimization algorithm, the present invention can generate regular quadrilateral meshes, reduce the calculation errors in seismic numerical simulation, and make the finite element analysis results more accurate and reliable.

[0046] (2) The present invention eliminates topological errors and deformed quadrilaterals such as concave quadrilaterals in the grid by optimizing the grid shape, thereby improving the stability of numerical calculations. The regular grid shape reduces numerical instability and oscillation phenomena during the numerical calculation process, improves the convergence of the algorithm, and ensures the stability of seismic numerical simulation.

[0047] (3) In the present invention, the side length variance of the edges around the grid nodes is incorporated into the loss function, so that the side length of the quadrilateral grid is evenly distributed and the edges that are too small in the grid are eliminated. In earthquake numerical simulation, time needs to be discretized. The smaller the grid side length, the smaller the discretization time interval, and the greater the required computational cost. The present invention optimizes the grid side length to make it more uniform, avoids the appearance of too small side length in the grid, and thus reduces the computational time and memory consumption.

[0048] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the embodiments of the present invention are specifically cited below and described in detail with reference to the attached drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0050] Figure 1 It is a specific implementation flow chart of the present invention;

[0051] Figure 2 is a data structure diagram for storing grid node information in an embodiment of the present invention;

[0052] Figure 3 is a data structure diagram for storing quadrilateral information in an embodiment of the present invention;

[0053] Figure 4 Schematic diagram of four types of nodes to be optimized in an embodiment of the present invention;

[0054] Figure 5 Schematic diagram of node solution space in an embodiment of the present invention;

[0055] Figure 6 This is a comparison diagram before and after optimization of some areas in an embodiment of the present invention. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all of the embodiments.

[0057] The embodiment of the present invention provides a method for optimizing the quadrilateral mesh of a geological model based on a particle swarm optimization algorithm. The actual data is the quadrilateral mesh data of a geological model collected in a certain area of ​​China. The mesh data consists of 97687 mesh nodes and 97058 quadrilaterals. In this example, the present invention effectively solves the problems of low calculation accuracy, low efficiency, poor stability, etc. of the geological model in finite element calculation due to the poorly shaped quadrilaterals.

[0058] Please refer to Figure 1 The specific optimization method of this embodiment includes the following steps:

[0059] S101, reading geological data including quadrilateral grids by parsing the divided geological model file, and storing the data in an internal data structure composed of a grid node table and a quadrilateral connection table.

[0060] In this embodiment, grid information is loaded from a standard grid file of a geological finite element model, the grid information including the numbers and coordinates of grid nodes, the numbers of quadrilaterals and the numbers of their constituent nodes, and the numbers of grid nodes on the boundary of the geological finite element model; all grid node information is traversed, and three types of information, namely, the numbers of grid nodes, the coordinates of grid nodes, and whether the grid nodes are on the model boundary, are stored in a grid node table; all quadrilateral information is traversed, and the numbers of the quadrilaterals and the numbers of the four constituent nodes are stored in the quadrilateral table, and the number of the quadrilateral is stored in the information of the four constituent nodes in the grid node table.

[0061] S102, traverse all nodes in the grid node table except the boundary nodes, add the nodes to be optimized to the node table to be optimized L, if the length of the node table to be optimized L is 0, jump to step S113, otherwise continue to execute step S103.

[0062] In this embodiment, it is necessary to detect whether the node needs to be optimized, that is, to detect whether the node is a node to be optimized.

[0063] There are four types of nodes to be optimized:

[0064] Category a: nodes that cause topological errors;

[0065] Category b: nodes that cause concave quadrilaterals;

[0066] Category C: nodes that cause quadrilaterals to degenerate into triangles;

[0067] Category D: nodes that cause angles that are too large or too small (the threshold for measuring whether the angle is too large or too small is set by the user. In this example, angles less than 22° are considered too small, and angles greater than 160° are considered too large).

[0068] The nodes in the node list L to be optimized will be used in the subsequent optimization algorithm, and the optimization order priority of the four types of nodes to be optimized decreases from a to d.

[0069] Methods for determining whether a node causes a topological error include:

[0070] Find all quadrilaterals with the node as a vertex, and find all line segments on the quadrilateral that are not connected to the node. All line segments can enclose a closed area. Whether the node generates a topological error depends on whether the node is inside the closed area.

[0071] Take this node as the starting point and draw a ray in the negative direction of the y axis;

[0072] Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node has a topological error.

[0073] Methods for determining whether a node forms a concave quadrilateral include:

[0074] Find all nodes connected to the node, connect these nodes to form a closed area, and detect whether the node will lead to the appearance of a concave quadrilateral, which depends on whether the node is inside the closed area;

[0075] Take this node as the starting point and draw a ray in the negative direction of the y axis;

[0076] Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node causes the appearance of a concave quadrilateral.

[0077] S103, initialize node number i=1.

[0078] S104, determine whether i is greater than the length of the node table L to be optimized, if so, jump to step S102, otherwise continue to execute step S105.

[0079] S105, select the i-th node N from the node list L to be optimized i , determine the node N i The solution space of i The initial positions of all particles in the particle swarm are initialized in the solution space (during the particle optimization process, the movement range of the particles is also carried out in the solution space of the node. In this embodiment, there are 500 particles in the particle swarm), and the loss value of the initial position is calculated according to the loss function of the optimization target, and then the individual optimal position and the global optimal position are initialized, and it is ensured that all particles will not leave the solution space during the movement.

[0080] In this embodiment, for node N i , find all nodes N i Connected nodes, these nodes are connected to form a closed area, that is, node N i The solution space is the range in which particles can be initialized and moved. If the particle movement position exceeds the solution space of the node, a quadrilateral topology error will occur or a concave quadrilateral will appear.

[0081] The method for initializing the individual optimal position and the global optimal position is: calculate the loss value of the initial position of each particle according to the loss function of the optimization target, and initialize the initial position with the lowest loss value as the global optimal position; for a certain particle, use the initial position of the particle as the initial individual optimal position of the particle.

[0082] Calculate the node N i Ways to lose value at a position include:

[0083] Step 1: Calculate the number of nodes N i is the loss value of all quadrilaterals at the vertex, and the formula is as follows:

[0084]

[0085] Among them, loss quad is the loss value of a single quadrilateral; TopologyError indicates that there is a topological error in the quadrilateral, and two of the four sides are crossed. TopologyErrorLoss is the loss value of the quadrilateral in this case; ConcaveError indicates that the quadrilateral is a concave quadrilateral, and ConcaveErrorLoss is the loss value of the quadrilateral in this case; AngleError indicates that the vertex angle of the quadrilateral is too large or too small, and AngleErrorLoss is the loss value of the quadrilateral in this case; angle i is the i-th interior angle of the quadrilateral; len i is the ith ordered edge of the quadrilateral, Represents the ratio of two opposite sides in a quadrilateral; the loss values ​​in the four categories in the loss value calculation formula decrease from top to bottom.

[0086] In this example, the value of TopologyErrorLoss is set to 999999.0, the value of ConcaveErrorLoss is set to 888888.0, and the value of AngleErrorLoss is set to 777777.0.

[0087] Step 2: Calculate the number of nodes N i The loss value (the surrounding nodes N i is the loss value of all quadrilaterals of the vertex, plus the surrounding nodes N i is the variance of the length of the edge at the endpoint), the formula is as follows:

[0088]

[0089] Among them, loss node is the loss value of the node, n is the number of quadrilaterals with the node as the vertex, is the loss value of the ith quadrilateral, m is the number of edges with this node as the endpoint, len j and len k are the lengths of the j-th and k-th edges.

[0090] In this embodiment, the method for determining whether two edges intersect is as follows:

[0091] Step 1: The two opposite sides of the quadrilateral are AB and CD. Calculate the projection of the two sides on the x-axis. If the projections do not overlap, it means that the two sides do not intersect. Otherwise, further judgment is made.

[0092] Step 2: Calculate the vector Calculating vectors and The cross product of can get the position of point C relative to line segment AB, and calculate the vector and The cross product of gives the position of point D relative to line segment AB.

[0093] Step 3: If the two cross product results have the same sign, it means that point C and point D are on the same side of line segment AB, and there is no intersection between AB and CD. If the two cross product results have different signs, it means that there is an intersection.

[0094] The method to determine whether a quadrilateral is a concave quadrilateral is as follows:

[0095] N1, N2, N3, N4 are the four ordered vertices of the quadrilateral. and is the area of ​​the triangle formed by the three vertices N1, N2, and N3, is the area of ​​the triangle formed by the three vertices N1, N4, and N3, is the area of ​​the triangle formed by the three vertices N2, N1, and N4, is the area of ​​the triangle formed by the three vertices N2, N3, and N4. If the two calculation results are equal, it means that the quadrilateral is a convex quadrilateral, otherwise it is a concave quadrilateral.

[0096] S106. Initialize iteration round j=1.

[0097] S107, determine whether the iteration round j is greater than the set maximum optimization iteration number M, if so, assign a value to i, i=i+1, and jump to step S104, otherwise continue to execute step S108.

[0098] S108, in the jth iteration, the particle number is initialized to k = 1, and the loss value of the optimal position of the individual particle is LocalLoss k .

[0099] S109, determine whether the particle number k is greater than the set total number of particles Q, if so, assign a value to j, j=j+1, and jump to step S107, otherwise continue to execute step S110.

[0100] S110, randomly generate particles P k The speed of the particle is calculated and the particle position is updated (to ensure that the particle position is confined to the solution space). During the position update process, the particle's approach to the individual optimal position and the global optimal position is balanced based on the weight parameter, and the particle P is calculated. k The loss value Loss at the latest position.

[0101] In this example, the weight of the particle's own velocity is 0.4, and the weights of approaching the individual optimal position and the global optimal position are both 0.3.

[0102] S111, if the loss value Loss at the latest position is less than the loss value LocalLoss at the optimal position of the particle individual k , then update the current position to the individual optimal position of particle k, and continue to execute step S112, otherwise the kth particle is optimized, assign a value to k, k=k+1, and return to step S109.

[0103] S112, if the loss value Loss at the latest position is less than the loss value at the global optimal position of the particle swarm, then update the current position to the global optimal position of the particle swarm, the kth particle is optimized, assign a value to k, k=k+1, and return to step S109.

[0104] S113, saving the optimized quadrilateral grid data of the geological model.

[0105] Figure 2 This is a data structure diagram for storing grid node information in this embodiment. The grid node information includes the index of the grid node, the X coordinate and Y coordinate of the node, and whether it is a node on the boundary. Q1, Q2, Q3, etc. are the indices of the quadrilateral where the node is located. Figure 3 The data structure diagram for storing the mesh quadrilateral in this embodiment includes the index of the quadrilateral, and N1, N2, N3, etc. are the indexes of the nodes constituting the quadrilateral. Figure 4 Schematic diagram of four types of nodes to be optimized in this embodiment, where Figure (a) shows a node that generates a topological error, Figure (b) shows a node whose inner angle of a quadrilateral exceeds 180°, Figure (c) shows a node in a quadrilateral that degenerates into a triangle due to an inner angle of 180°, and Figure (d) shows a node whose inner angle of a quadrilateral does not exceed 180° but exceeds the maximum inner angle initially set. Figure 5 Schematic diagram of the grid node solution space in this embodiment. The closed area surrounded by the dotted line is the solution space of the hollow node in the figure. If the node position exceeds the solution space range, topological errors or concave quadrilaterals will occur. Figure 6 These are comparisons of some areas of this example before and after optimization. The large image above is the quadrilateral mesh image before optimization, and the large image below is the quadrilateral mesh image after optimization. The small images above areas A, B, C, and D are areas that need to be optimized in the quadrilateral mesh before optimization, which respectively contain a, b, c, and d nodes in the nodes to be optimized. The small images below areas A, B, C, and D are areas after the corresponding node positions in the small images above are optimized.

[0106] This embodiment can illustrate the geological model quadrilateral mesh optimization method based on the particle swarm optimization algorithm described in the present invention, which effectively solves the problems of low calculation accuracy, low efficiency, poor stability, etc. caused by poorly shaped quadrilaterals in actual seismic finite element calculations, and at the same time provides an efficient and general optimization method for similar quadrilateral mesh optimization in the future.

[0107] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for optimizing quadrilateral meshes of geological models based on particle swarm optimization algorithm, characterized in that: The following steps are involved: S101, read the geological model grid file obtained by the dissection, and store the read file data in an internal data structure composed of a grid node table and a quadrilateral table; S102, traverse all nodes in the grid node table except the boundary nodes, add the nodes to be optimized to the node table to be optimized L, if the length of the node table to be optimized L is 0, jump to step S113, otherwise continue to execute step S103; S103, start traversing the node table L to be optimized, and initialize the node number i=1; S104, determine whether i is greater than the length of the node table L to be optimized, if so, jump to step S102 and perform the next round of optimization on the entire grid, otherwise continue to step S105; S105, select the i-th node N from the node list L to be optimized i , determine the node N i The solution space of the particle swarm is initialized in the solution space, and the loss value of the initial position is calculated according to the loss function of the optimization target. Then, the individual optimal position and the global optimal position are initialized, and it is ensured that all particles will not leave the solution space during the movement. S106, start iterative optimization of node N i Position, initialize the iteration round j = 1; S107, determine whether the iteration round j is greater than the set maximum optimization iteration number M, if so, assign a value to i, i=i+1, and jump to step S104, otherwise continue to execute step S108; S108, in the jth iteration, the particle number is initialized to k = 1, and the loss value of the optimal position of the individual particle is LocalLoss k ; S109, determine whether the particle number k is greater than the set total number of particles Q, if so, assign a value to j, j=j+1, and jump to step S107, otherwise continue to execute step S110; S110, randomly generate particles P k The speed of the particle is calculated and the position of the particle is updated. During the position update process, the particle's approach to the individual optimal position and the global optimal position is balanced based on the weight parameter, and the particle P is calculated. k The loss value Loss at the latest position; S111, if the loss value Loss at the latest position is less than the loss value LocalLoss at the optimal position of the particle individual k , then update the current position to the individual optimal position of particle k, and continue to execute step S112, otherwise the optimization of the kth particle is completed, assign a value to k, k=k+1, and return to step S109; S112, if the loss value Loss at the latest position is less than the loss value at the global optimal position of the particle swarm, then update the current position to the global optimal position of the particle swarm, the kth particle is optimized, assign a value to k, k=k+1, and return to step S109; S113, saving the optimized quadrilateral grid data of the geological model.

2. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 1, characterized in that: Step S101 specifically includes: Parse and load mesh information from the standard mesh file used by the seismic finite element model. The mesh information includes the numbers and coordinates of mesh nodes, the numbers of quadrilaterals and their constituent nodes, and the numbers of mesh nodes on the boundaries of the seismic finite element model. Traverse all grid node information, and store three types of information, namely, the number of the grid node, the coordinates of the grid node, and whether the grid node is at the boundary of the seismic finite element model, in the grid node table; All quadrilateral information is traversed, and the quadrilateral number and the numbers of the four constituent nodes are stored in the quadrilateral table, and the quadrilateral number is stored in the information of the four constituent nodes in the grid node table.

3. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 1, characterized in that: In step S102, the nodes to be optimized include four categories, a, b, c and d. The four categories of nodes are nodes that cause topological errors, nodes that cause concave quadrilaterals, nodes that cause quadrilaterals to degenerate into triangles, and nodes that cause angles that are too large or too small. The optimization order priority of the four categories of nodes to be optimized decreases from a to d.

4. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 3, characterized in that: Methods for determining whether a node causes a topological error include: Find all quadrilaterals with the node as a vertex, and find all line segments on the quadrilateral that are not connected to the node. All line segments can enclose a closed area. Take this node as the starting point and draw a ray in the negative direction of the y axis; Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node has a topological error.

5. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 3, characterized in that: Methods for determining whether a node forms a concave quadrilateral include: Find all nodes connected to this node, and connect these nodes to form a closed area; Take this node as the starting point and draw a ray in the negative direction of the y axis; Find the number of intersections between the ray and the boundary of the closed area. If the number of intersections is an even number, it means that the node causes the appearance of a concave quadrilateral.

6. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 3, characterized in that: Calculate the node N i Ways to lose value at a position include: Step 1: Calculate the number of nodes N i is the loss value of all quadrilaterals at the vertex, and the formula is as follows: Among them, loss quad is the loss value of a single quadrilateral; TopologyError indicates that there is a topological error in the quadrilateral, and two of the four sides are crossed. TopologyErrorLoss is the loss value of the quadrilateral in this case; ConcaveError indicates that the quadrilateral is a concave quadrilateral, and ConcaveErrorLoss is the loss value of the quadrilateral in this case; AngleError indicates that the vertex angle of the quadrilateral is too large or too small, and AngleErrorLoss is the loss value of the quadrilateral in this case; angle i is the i-th interior angle of the quadrilateral; len i is the ith ordered edge of the quadrilateral, It represents the ratio of two opposite sides in a quadrilateral; Step 2: Calculate the number of nodes N i The loss value is as follows: Among them, loss node is the loss value of the node, n is the number of quadrilaterals with the node as the vertex, is the loss value of the ith quadrilateral, m is the number of edges with this node as the endpoint, len j and len k are the lengths of the j-th and k-th edges.

7. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 1, characterized in that: In step S105, for node N i , find all nodes N i Connected nodes, these nodes are connected to form a closed area, that is, node N i solution space; calculate the loss value of the initial position of each particle according to the loss function of the optimization objective, and initialize the initial position with the lowest loss value as the global optimal position; for a certain particle, take the initial position of the particle as the initial individual optimal position of the particle.

8. The method for optimizing quadrilateral mesh of geological model based on particle swarm optimization algorithm according to claim 1, characterized in that: For a particle p k , the optimization process includes: Calculate particle p k The loss value Loss at the latest position; If the loss value Loss is greater than or equal to its individual optimal position, then particle p k Complete the optimization, otherwise update the current position to particle p k The individual optimal position of If the loss value Loss is greater than or equal to the global optimal position, then particle p k Complete the optimization, otherwise update the current position to the global optimal position and complete the particle p k Optimization.

Images