Finite difference time domain modeling method and system for undulating terrain and irregular anomaly

Abstract

The present invention belongs to the technical field of transient electromagnetic three-dimensional forward modeling. It provides a finite difference time domain modeling method and system for undulating terrain and irregular anomaly bodies, models complex target body triangular meshes, uses improved parallel ray tracing technology to generate intersection information of triangular meshes and FDTD meshes, and uses conformal mesh technology to calculate the equivalent conductivity of edges intersecting with complex interfaces, effectively overcoming the step approximation error caused by the existing FDTD method and achieving fine characterization of complex geological structures.

Description

Technical Field

[0001] The invention relates to the technical field of transient electromagnetic three-dimensional forward modeling, and in particular to a finite difference time-domain modeling method and system for undulating terrain and irregular anomaly bodies. Background Art

[0002] The statements in this section merely provide background art related to the present invention and do not necessarily constitute prior art.

[0003] Finite-Difference Time-Domain (FDTD) is one of the commonly used methods for transient electromagnetic three-dimensional forward modeling. This method can accurately reflect the response law of three-dimensional geological structures under transient electromagnetic field excitation. It has the advantages of simple programming, low memory consumption and high computational efficiency.

[0004] The actual application scenarios of transient electromagnetic detection are complex, with rugged terrain and irregular shapes of detection targets. When using the conventional FDTD method to perform three-dimensional forward simulation calculations on complex geological structures, it is necessary to use structured grids to discretize complex terrain or irregular anomalies, which results in a step approximation error, and this error increases with the increase of grid size; using a smaller grid division model can effectively reduce the step approximation error, but it also means that the number of grids after the calculation model is discretized increases, which inevitably leads to an increase in memory usage and an extension of calculation time. Summary of the invention

[0005] In order to address the deficiencies of the prior art, the present invention provides a finite difference time-domain modeling method and system for undulating terrain and irregular anomalies, models the triangular mesh of complex target bodies, uses improved parallel ray tracing technology to generate the intersection information of the triangular mesh and the FDTD mesh, and uses conformal mesh technology to calculate the equivalent conductivity of the edges intersecting with the complex interface, effectively overcoming the step approximation error caused by the existing FDTD method and achieving a fine characterization of complex geological structures.

[0006] In order to achieve the above object, the present invention adopts the following technical solution:

[0007] In a first aspect, the present invention provides a finite difference time domain modeling method for undulating terrain and irregular anomaly bodies.

[0008] A finite difference time domain modeling method for undulating terrain and irregular anomaly bodies includes the following processes:

[0009] Generate a surface model based on the point cloud data of undulating terrain and irregular anomalies, and use triangular meshes to divide the surface model to obtain multiple triangular face elements;

[0010] Modify the normal direction of each triangular face element to obtain the triangular face element position information of the undulating terrain and the irregular anomaly body respectively;

[0011] Import the triangular surface element position information into the FDTD three-dimensional forward model and calculate the normal vector of each triangular surface element;

[0012] The FDTD grid is used to divide the calculation model. The FDTD grid is regarded as a ray that shoots from one side of the model to the other side along the three coordinate axes of x, y and z. The starting point and direction vector of each ray are recorded. The minimum storage ray-triangle intersection algorithm is used to obtain the intersection point of each triangle face element and the ray.

[0013] Determine the entry and exit attributes of all intersection points based on the inner product of the normal vector of each triangle face and the ray;

[0014] According to the in and out properties of all intersection points, the conformal grid is used to calculate the equivalent conductivity at the edges of the conformal grid;

[0015] The equivalent conductivity is substituted into the electric field iteration equation of the FDTD algorithm, and transient electromagnetic three-dimensional forward calculation is performed to obtain the FDTD modeling results.

[0016] As a further limitation of the first aspect of the present invention, before generating a surface model based on the point cloud data of the undulating terrain and the irregular anomaly, the following process is included:

[0017] Extract the latitude, longitude and elevation data of the terrain in the area to be measured from satellite remote sensing data and convert them into an engineering coordinate system;

[0018] The surface model is constructed in the engineering coordinate system.

[0019] As a further limitation of the first aspect of the present invention, modifying the normal direction of the triangular mesh includes: ensuring that the normal directions of the triangular face elements after the terrain surface is divided uniformly point to the air or the stratum, and the normal directions of the triangular face elements after the anomaly surface is divided uniformly point to the inside or outside of the anomaly.

[0020] As a further limitation of the first aspect of the present invention, the normal vector of each triangular face element is calculated ,include:

[0021] Assume that the three vertices of the triangle are , and :

[0022] ;

[0023] .

[0024] As a further limitation of the first aspect of the present invention, the minimum stored ray-triangle intersection algorithm includes:

[0025] Assume the ray equation is: , where O is the starting point of the ray and t is the extension length of the ray. D is the direction vector of the ray, along the x, y, and z axes are (1,0,0), (0,1,0), and (0,0,1), respectively, with a triangular face element , the coordinates of any point in the surface element can be expressed using the barycentric coordinate system as follows: ,in, , and Triangular face The coordinates of the three vertices;

[0026] Combining the ray and the coordinate equation of any point on the triangular surface element, we get: ;

[0027] Modify it to matrix multiplication format and get: ;

[0028] Cramer's rule can be solved to obtain: ;

[0029] like t, u, v Requirements: t >0, 0≤ u ≤1, 0≤ v ≤1 and , then the current triangle element intersects with the FDTD grid, and t is substituted into the ray equation to obtain the Cartesian coordinates of the intersection point. After the calculation is completed, the intersection points in the x, y and z directions are stored in three arrays respectively.

[0030] As a further limitation of the first aspect of the present invention, determining the entry and exit attributes of all intersection points according to the inner product of the normal vector of each triangular face element and the ray includes:

[0031] The direction vector of the x-ray is =(1, 0, 0), the inner product of the x-direction ray and the normal vector passing through the triangle surface element. If the inner product is less than 0, the projection of the triangle surface element normal vector on the ray direction vector is negative, then the intersection attribute of the triangle surface element and the x-direction ray is the entry point; if the inner product is greater than 0, the projection of the triangle surface element normal vector on the ray direction vector is positive, then the intersection attribute of the triangle surface element and the x-direction ray is the exit point.

[0032] As a further limitation of the first aspect of the present invention, for undulating terrain, the part between the penetration point and the exit point is the anomaly, and the part between the exit point and the next penetration point is air; for irregular anomalies, the part between the penetration point and the exit point is the anomaly, and the rest is the stratum.

[0033] As a further limitation of the first aspect of the present invention, when a ray passes through a common edge or common vertex of a triangular face, the entry and exit attributes of each triangular face intersection where the common edge or common vertex is located are determined. If the entry and exit attributes of all intersections are consistent, it means that the ray enters or exits an irregular anomaly or undulating terrain, and only one intersection needs to be retained. If the entry and exit attributes of the intersection are inconsistent, it means that the ray enters or exits an anomaly or terrain, and the medium on both sides of the intersection has not changed, and the intersection is not retained.

[0034] As a further limitation of the first aspect of the present invention, using a conformal grid to calculate the equivalent conductivity at the edges of the conformal grid includes:

[0035] The target body surface after triangular facet segmentation divides the Yee cell into two spaces, which are filled with two media with different electrical parameters respectively. The electric field value is calculated at the center of the cell edge, and the equivalent conductivity is obtained by the weighted average calculation of the length occupied by different media on the corresponding edge.

[0036] In a second aspect, the present invention provides a finite difference time domain modeling system for undulating terrain and irregular anomalies.

[0037] A finite difference time domain modeling system for undulating terrain and irregular anomaly bodies, comprising:

[0038] The triangular face element generation unit is configured to: generate a surface model according to the point cloud data of the undulating terrain and the irregular anomaly, and divide the surface model using a triangular mesh to obtain a plurality of triangular face elements;

[0039] The position information generating unit is configured to: modify the normal direction of each triangular face element to obtain the triangular face element position information of the undulating terrain and the irregular abnormal body respectively;

[0040] The normal vector calculation unit is configured to: import the position information of the triangular surface element into the FDTD three-dimensional forward model, and calculate the normal vector of each triangular surface element;

[0041] The ray tracing unit is configured to: use FDTD grid to divide the calculation model, the FDTD grid is regarded as a ray shooting from one side of the model to the other side along the three coordinate axes of x, y and z, record the starting point and direction vector of each ray, and use the minimum storage ray and triangle intersection algorithm to obtain the intersection point of each triangle face element and the ray;

[0042] The entry and exit attribute judgment unit is configured to: determine the entry and exit attributes of all intersection points according to the inner product of the normal vector of each triangular face element and the ray;

[0043] The equivalent conductivity calculation unit is configured to: calculate the equivalent conductivity at the edges of the conformal grid using the conformal grid according to the in and out attributes of all intersection points;

[0044] The modeling result generating unit is configured to: substitute the equivalent conductivity into the electric field iteration equation of the FDTD algorithm, perform transient electromagnetic three-dimensional forward calculation, and obtain the FDTD modeling result.

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

[0046] The present invention realizes the refined modeling and meshing of terrain and anomalies, adopts the improved parallel ray tracing technology, realizes the fast intersection of triangular mesh and FDTD mesh; adopts the conformal mesh technology, effectively reduces the step approximation error brought by the traditional FDTD method in the numerical simulation of complex interfaces, and improves the numerical simulation accuracy of complex terrain and irregular anomalies in transient electromagnetic three-dimensional forward simulation.

[0047] Advantages of additional aspects of the present invention will be given in part in the following description, and in part will become obvious from the following description, or will be learned through practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] The accompanying drawings in the specification, which constitute a part of the present invention, are used to provide a further understanding of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute improper limitations on the present invention.

[0049] Figure 1 A schematic diagram of the latitude, longitude and elevation data of the terrain of the area to be measured provided in Example 1 of the present invention;

[0050] Figure 2 A schematic diagram of triangulated meshing of undulating terrain provided in Example 1 of the present invention;

[0051] Figure 3 Vertex coordinate information of the .dat file provided in Example 1 of the present invention;

[0052] Figure 4 The metadata information of the .dat file provided in Example 1 of the present invention;

[0053] Figure 5 A schematic diagram showing that the order of vertices and the direction of normals of triangular face elements provided in Example 1 of the present invention conform to the right-hand screw rule;

[0054] Figure 6 A schematic diagram of the intersection of a ray and a triangular face provided in Embodiment 1 of the present invention;

[0055] Figure 7 A schematic diagram of the intersection point provided in Embodiment 1 of the present invention being located at the common side of a triangle and the ray passing through an abnormal body;

[0056] Figure 8A schematic diagram of the embodiment 1 of the present invention in which the intersection point is located on the common side of the triangle but the ray does not pass through the abnormal body;

[0057] Fig. 9 A schematic diagram of equivalent medium parameters in a conformal grid of a Yee cell provided in Example 1 of the present invention;

[0058] Fig.10 Provided in Example 1 of the present invention xoz Schematic diagram of calculation of equivalent conductivity of surface;

[0059] Fig.11 A schematic diagram of a spherical anomaly model in a uniform half space provided in Example 1 of the present invention;

[0060] Fig.12 Three methods provided for Example 1 of the present invention Response comparison chart;

[0061] Fig.13 A schematic diagram of a finite difference time domain modeling system for undulating terrain and irregular anomaly bodies provided in Example 2 of the present invention. DETAILED DESCRIPTION

[0062] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0063] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0064] In the absence of conflict, the embodiments of the present invention and the features of the embodiments may be combined with each other.

[0065] Embodiment 1:

[0066] The dielectric conformal grid technology uses the electrical parameter mapping between two sets of grids to perform special processing on the surface grid of complex structures. It can reduce the step approximation error without refining the grid, thereby achieving a fine characterization of complex structures. At present, in the field of transient electromagnetic three-dimensional forward modeling, conformal grid technology is rarely used, and the boundary characterization of complex terrain and irregular anomalies still relies heavily on step approximation.

[0067] In order to realize conformal mesh technology, it is necessary to first fine-tune the modeling of complex structures. The triangular facet modeling method is widely used in electromagnetic modeling because it can better approximate the actual situation of complex terrain and irregular anomaly surfaces, and the model file information has the advantages of simple and intuitive, easy information acquisition, etc. At present, the modeling technology of commercial software is very mature. Among them, the GiD (Group Identification) platform is a general, adaptable, user-friendly and powerful pre- and post-processing system that can be used for geometric simulation, data input, mesh generation, core algorithm iteration and visualization of calculation results. GiD can quickly establish geometric models or import model files in IGES, DXF, OBJ, ACIS and other formats to realize structured or unstructured mesh generation of models. It can provide three-dimensional visualized geometric modeling and Yee cell format mesh generation for FDTD numerical simulation, which provides a platform support for modeling research.

[0068] The present invention models the triangular mesh of complex target bodies through GiD, uses improved parallel ray tracing technology to generate the intersection information of the triangular mesh and the FDTD mesh, and adopts conformal mesh technology to calculate the equivalent conductivity of the edges intersecting with the complex interface, which effectively overcomes the step approximation error caused by the existing FDTD method and realizes the fine characterization of complex geological structures.

[0069] Specifically, this implementation provides a finite difference time domain modeling method for undulating terrain and irregular anomaly bodies, including the following process:

[0070] S1: Refined modeling of terrain and irregular anomalies. Extract the latitude, longitude and elevation data of the survey area terrain from satellite remote sensing data and convert them into engineering coordinate system; point cloud data of irregular anomalies such as real ore bodies and underground caves can be generated by existing drilling geological data, laser radar scanning and other means; the accuracy of terrain and anomaly characterization is related to the number of sampling points. The denser the sampling points, the finer the terrain characterization and the larger the data volume.

[0071] S2: Import the point cloud data describing the terrain and anomalies into GiD software, generate the surface model of the structure, and use triangulated meshing to divide the surface model, such as Figure 1 and Figure 2 As shown in the figure, modify the normal direction of the triangular mesh to ensure that the normal direction of the triangular face element after the terrain surface is uniformly pointed to the air or the stratum, and the normal direction of the triangular face element after the abnormal body surface is uniformly pointed to the inside or outside of the abnormal body. Output the .dat files storing the position information of the undulating terrain and the irregular abnormal body triangular mesh respectively, such as Figure 3 and Figure 4 As shown, Figure 3 The number of triangle face fixed points, the number of triangle face elements, vertex numbers and vertex coordinates are given; Figure 4The triangular face elements and corresponding vertices are given, including the face element number, vertex 1, vertex 2 and vertex 3.

[0072] The order of the three vertices of the output triangle element and the direction of the normal line conform to the right-hand screw rule. If the normal line is uniformly directed outward, the vertex numbering of the triangle element output by GiD is counterclockwise, such as Figure 5 As shown, specifying the normal direction of the triangular face element determines the output order of the three vertices of the face element, which gives an example of the face element number and the corresponding vertex 1, vertex 2 and vertex 3.

[0073] S3: Import the .dat file into the FDTD conformal mesh program, which will automatically read the triangle element and vertex information. Then calculate the normal vector of each triangle element. . Let the three vertices of the triangle be , whose normal vector can be calculated by vector cross product:

[0074] (1);

[0075] (2);

[0076] Move the four fingers of your right hand from Steering , the direction the thumb points to is the normal vector direction.

[0077] S4: Use uniform or non-uniform FDTD grid to divide the calculation model, where the FDTD grid can be regarded as a ray shooting from one side of the model to the other side along the three coordinate axes of x, y, and z. Record the starting point and direction vector of each ray, and use the improved parallel ray tracing technology to calculate the intersection coordinate information of the triangular face elements generated after the anomaly and terrain division and the ray.

[0078] The traditional method of determining whether a ray and a triangle intersect is divided into two steps: first, determining whether the ray and the plane where the triangle is located have an intersection point and the intersection coordinates, and then determining whether the intersection point is inside the triangle. This method has redundant steps and low computational efficiency. The present invention adopts a ray and triangle intersection algorithm with minimum storage, introduces parallel programming technology, and realizes a one-time batch solution of intersection points. The algorithm implementation steps are as follows:

[0079] Let the ray equation for:

[0080] (3);

[0081] in, O is the starting point of the ray, t is the extension length of the ray, Dis the direction vector of the ray, which is (1,0,0), (0,1,0) and (0,0,1) along the x, y and z axes respectively; there is a triangular face element , the coordinates of any point in the surface element can be expressed using the barycentric coordinate system as follows:

[0082] (4);

[0083] in, , and Triangular face The coordinates of the three vertices.

[0084] Combining the ray and the coordinate equation of any point on the triangular surface element, we get:

[0085] (5);

[0086] Rearranging formula (5) into matrix multiplication format yields:

[0087] (6);

[0088] Cramer's rule can be solved to obtain:

[0089] (7);

[0090] like t, u, v At the same time, the following conditions must be met: (1) ; (2) , ; (4) , indicating that the triangular face element intersects with the FDTD grid. The Cartesian coordinates of the intersection point can be obtained by substituting t into equation (3). After the calculation is completed, the intersection points in the x, y, and z directions are stored in three arrays respectively.

[0091] The intersection process of each triangle face element and the ray is independent and does not interfere with each other, which meets the conditions for parallel programming. Therefore, the CPU parallel programming model OpenMP is introduced to realize the parallelization of the ray tracing algorithm, and the calculation speed is significantly improved.

[0092] S5: When the ray passes through a complex terrain surface, there is generally only one intersection point along the z-axis. Along the z-axis, there is an air layer above the terrain interface and a stratum below the interface. The equivalent conductivity needs to be calculated at the edge of the conformal grid where the terrain interface is located.

[0093] When the ray passes through an abnormal body or through undulating terrain along the x- and y-axis directions, multiple pairs of intersection points will generally appear. Figure 6 Taking the ray in the x direction as an example, its direction vector is =(1, 0, 0), the three vertices of the penetrated triangle are A1, A2 and A3, the penetration point is P1, and the normal vector is ; The three vertices of the triangle element are B1, B2 and B3, the entry point is Q1, and the normal vector is .

[0094] Calculate the inner product with the normal vector of the triangle face element. If the inner product is less than 0, that is, the projection of the normal vector of the triangle face element on the ray direction vector is negative, then the intersection attribute of the triangle face element and the x-direction ray is the entry point; if the inner product is greater than 0, that is, the projection of the normal vector of the triangle face element on the ray direction vector is positive, then the intersection attribute of the triangle face element and the x-direction ray is the exit point. For undulating terrain, the part between the entry point and the exit point is the anomaly, and the part between the exit point and the next entry point is air; for an anomaly, the part between the entry point and the exit point is the anomaly, and the rest is the stratum.

[0095] S6: Determination of special intersections, including:

[0096] When the ray passes through the common edge or common vertex of the triangular face, the FDTD mesh may be tangent to the abnormal body or terrain. In this case, it is necessary to determine the entry and exit properties of each triangular face intersection where the common edge or common vertex is located. Figure 7 In the intersection case shown, if the entry and exit attributes of all intersection points are consistent, it means that the ray enters (exits) the abnormal body or terrain, and only one intersection point needs to be retained. Figure 7 In the figure, the four vertices of the triangle are D1, D2, D3 and D4, and the ray along the z direction is (0,0,1). If the in and out attributes of the intersection are inconsistent, it means that the medium on both sides of the intersection has not changed after entering or exiting the abnormal body or terrain. Figure 8 The intersection case shown is not considered, where the four vertices of the triangle are E1, E2, E3 and E4, and the ray along the z direction is (0,0,1).

[0097] S7: After determining the entry and exit attributes of all intersections, the conformal mesh technique is used to calculate the equivalent conductivity of the edges where the triangle mesh and the FDTD mesh intersect. Fig. 9 As shown in the figure, the target surface after triangular facet segmentation divides the Yee cell into two spaces, which are filled with two media with different electrical parameters. Since the calculation of the electric field value is located at the center of the cell edge, the equivalent conductivity is obtained by the weighted average calculation of the length occupied by different media on the corresponding edge.

[0098] In this implementation, Yee cell is a numerical method for simulating electromagnetic fields and wave phenomena. It is based on Maxwell equations and finite difference method. This method divides the space into small cubic units (cells). The field quantity in each unit is calculated in a discrete manner. The core idea of ​​Yee cell is to use finite difference to approximate partial differential equations, so as to solve the behavior of electromagnetic fields and wave phenomena.

[0099] In this implementation, preferably, the conductivity of the earth and the air is assigned using the undulating terrain as the interface, and the conductivity inside the anomaly is assigned using the surface of the anomaly as the interface.

[0100] consider xoz Surface conditions, such as Fig.10 As shown, , , and , , They represent the three-directional components of the electric field and magnetic field in the Yee cell space discretization. The edge medium 1 on the right side of the grid (with a conductivity of ) and medium 2 (conductivity is ) occupy lengths respectively and The lengths occupied by the lower edge medium 1 and medium 2 are and The edges where points A and B are located have no intersection with the triangular surface element, so there is no need to solve the equivalent conductivity, and the conductivity of medium 1 can be assigned; the equivalent conductivity of points C and D is calculated by weighted average of the corresponding lengths of the edges occupied by the two media. Fig.10 The equivalent conductivity of each edge is:

[0101] (8);

[0102] (9);

[0103] (10);

[0104] (11);

[0105] in, is the equivalent conductivity represented by edge A, is the equivalent conductivity represented by edge B, is the equivalent conductivity represented by edge C, is the equivalent conductivity represented by the edge D; the equivalent conductivity obtained in the above formula is substituted into the electric field iteration equation of the FDTD algorithm, and the magnetic field iteration equation does not need to be modified.

[0106] by Fig.11 The spherical anomaly model in the uniform half space shown in the figure is used as an example to verify the effectiveness of this modeling method. The size of the entire model is 2510m×2510m×2510m, and a uniform grid is used for partitioning, with a grid size of 10m ( = = ). The electrical conductivity of air, earth and sphere are S / m, 0.01S / m and 1S / m. The side length of the loop source is 90m, the transmission waveform adopts a trapezoidal wave processed by the switching function, the transmission current is 1A, and the receiving point is located at the center of the loop source. A three-dimensional sphere model was created using GiD with a radius of 50m. The sphere is divided into 1496 triangular face elements and 750 vertices. The sphere model is set directly below the center of the transmission source and buried at a depth of 25m.

[0107] The modeling method and the step approximation method proposed in the present invention are used to calculate the grid electrical parameters, and the FDTD method is used to calculate the above model and draw the receiving point. Response curve ( for z The finite element method uses tetrahedral meshes to simulate underground anomalies, so this method can more accurately divide the anomaly model than the FDTD method. Therefore, the numerical simulation of the model calculated by the finite element method is used as the benchmark solution. Fig.12 It can be seen that for the three methods of finite element method, time-domain finite difference method-step approximation, time-domain finite difference method-conformal grid, the numerical simulation results after processing with conformal grid technology are in good agreement with the benchmark solution. The traditional step approximation method shows a large deviation from the benchmark solution in the late observation time, which verifies the effectiveness of the modeling method proposed in this invention.

[0108] Embodiment 2:

[0109] like Fig.13 As shown, this implementation provides a finite difference time domain modeling system for undulating terrain and irregular anomaly bodies, including:

[0110] The triangular face element generation unit is configured to: generate a surface model according to the point cloud data of the undulating terrain and the irregular anomaly, and divide the surface model using a triangular mesh to obtain a plurality of triangular face elements;

[0111] The position information generating unit is configured to: modify the normal direction of each triangular face element to obtain the triangular face element position information of the undulating terrain and the irregular abnormal body respectively;

[0112] The normal vector calculation unit is configured to: import the position information of the triangular surface element into the FDTD three-dimensional forward model, and calculate the normal vector of each triangular surface element;

[0113] The ray tracing unit is configured to: use FDTD grid to divide the calculation model, the FDTD grid is regarded as a ray shooting from one side of the model to the other side along the three coordinate axes of x, y and z, record the starting point and direction vector of each ray, and use the minimum storage ray and triangle intersection algorithm to obtain the intersection point of each triangle face element and the ray;

[0114] The entry and exit attribute judgment unit is configured to: determine the entry and exit attributes of all intersection points according to the inner product of the normal vector of each triangular face element and the ray;

[0115] The equivalent conductivity calculation unit is configured to: calculate the equivalent conductivity at the edges of the conformal grid using the conformal grid according to the in-out properties of all intersections;

[0116] The modeling result generating unit is configured to: substitute the equivalent conductivity into the electric field iteration equation of the FDTD algorithm, perform transient electromagnetic three-dimensional forward calculation, and obtain the FDTD modeling result.

[0117] The specific working process of each of the above units is shown in the method in Example 1, which will not be repeated here.

[0118] It is understandable that the above-mentioned units can be separately or completely combined into one or several other units to constitute, or one (some) of the units can be further divided into multiple smaller units in function to constitute, which can achieve the same operation without affecting the realization of the technical effects of the embodiments of the present application. The above-mentioned units are divided based on logical functions. In practical applications, the functions of one unit can also be implemented by multiple units, or the functions of multiple units can be implemented by one unit. In other embodiments of the present application, the system may also include other units. In practical applications, these functions can also be implemented with the assistance of other units, and can be implemented by the collaboration of multiple units.

[0119] According to another embodiment of the present application, the system described in this embodiment can be constructed, and the method of Example 1 of the present application can be implemented by running a computer program (including program code) capable of executing the steps involved in the corresponding method described in Example 1 on a general-purpose computing device such as a computer, which includes processing elements and storage elements such as a central processing unit (CPU), random access memory (RAM), and read-only memory (ROM). The computer program can be recorded on, for example, a computer-readable recording medium, and loaded into the above-mentioned computing device through the computer-readable recording medium and run therein.

[0120] 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 finite difference time domain modeling method for undulating terrain and irregular anomaly bodies, characterized in that: The process includes: Generate a surface model based on the point cloud data of undulating terrain and irregular anomalies, and use triangular meshes to divide the surface model to obtain multiple triangular face elements; Modify the normal direction of each triangular face element to obtain the triangular face element position information of the undulating terrain and the irregular anomaly body respectively; Import the triangular surface element position information into the FDTD three-dimensional forward model and calculate the normal vector of each triangular surface element; The FDTD grid is used to divide the calculation model. The FDTD grid is regarded as a ray that shoots from one side of the model to the other side along the three coordinate axes of x, y and z. The starting point and direction vector of each ray are recorded. The minimum storage ray-triangle intersection algorithm is used to obtain the intersection point of each triangle face element and the ray. Determine the entry and exit attributes of all intersection points based on the inner product of the normal vector of each triangle face and the ray; According to the in and out properties of all intersection points, the conformal grid is used to calculate the equivalent conductivity at the edges of the conformal grid; Substitute the equivalent conductivity into the electric field iteration equation of the FDTD algorithm, perform transient electromagnetic three-dimensional forward calculation, and obtain the FDTD modeling results; According to the inner product of the normal vector of each triangle face element and the ray, the entry and exit attributes of all intersection points are determined, including: The direction vector of the x-ray is =(1, 0, 0), the inner product of the x-direction ray and the normal vector passing through the triangle surface element. If the inner product is less than 0, the projection of the triangle surface element normal vector on the ray direction vector is negative, then the intersection attribute of the triangle surface element and the x-direction ray is the entry point; if the inner product is greater than 0, the projection of the triangle surface element normal vector on the ray direction vector is positive, then the intersection attribute of the triangle surface element and the x-direction ray is the exit point.

2. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 1, characterized in that: Before generating a surface model based on point cloud data of undulating terrain and irregular anomalies, the following processes are included: Extract the latitude, longitude and elevation data of the terrain in the area to be measured from satellite remote sensing data and convert them into an engineering coordinate system; The surface model is constructed in the engineering coordinate system.

3. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 1, characterized in that: Modify the normal direction of the triangular mesh, including: ensuring that the normal directions of the triangular face elements after the terrain surface is divided uniformly point to the air or the stratum, and the normal directions of the triangular face elements after the abnormal body surface is divided uniformly point to the inside or outside of the abnormal body.

4. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 1, characterized in that: Calculate the normal vector for each triangle face ,include: Assume that the three vertices of the triangle are , and : ； 。 5. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 1, characterized in that: Minimal-memory ray-triangle intersection algorithm, including: Assume the ray equation is: , where O is the starting point of the ray and t is the extension length of the ray. D is the direction vector of the ray, along the x, y, and z axes are (1,0,0), (0,1,0), and (0,0,1), respectively, with a triangular face element , the coordinates of any point in the surface element can be expressed using the barycentric coordinate system as follows: ,in, , and Triangular face The coordinates of the three vertices; Combining the ray and the coordinate equation of any point on the triangular surface element, we get: ; Modify it to matrix multiplication format and get: ; Cramer's rule can be solved to obtain: ; like t, u, v Requirements: t >0, 0≤ u ≤1, 0≤ v ≤1 and , then the current triangle element intersects with the FDTD grid, Substitute the ray equation to get the Cartesian coordinates of the intersection point. After the calculation is completed, the intersection points in the x, y and z directions are stored in three arrays respectively.

6. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 1, characterized in that: For undulating terrain, the part between the penetration point and the exit point is the anomaly, and the part between the exit point and the next penetration point is air; for irregular anomaly, the part between the penetration point and the exit point is the anomaly, and the rest is the stratum.

7. The finite difference time domain modeling method for undulating terrain and irregular anomaly according to claim 6, characterized in that: When a ray passes through a common edge or common vertex of a triangular face, the entry and exit properties of each triangular face intersection where the common edge or common vertex is located are determined. If the entry and exit properties of all intersections are consistent, it means that the ray enters or exits an irregular anomaly or undulating terrain, and only one intersection needs to be retained. If the entry and exit properties of the intersection are inconsistent, it means that the ray enters or exits an anomaly or terrain, and the medium on both sides of the intersection has not changed, so the intersection is not retained.

8. The finite difference time domain modeling method for undulating terrain and irregular anomaly bodies according to claim 1, characterized in that: Use conformal mesh to calculate the equivalent conductivity at the edges of the conformal mesh, including: The target body surface after triangular facet segmentation divides the Yee cell into two spaces, which are filled with two media with different electrical parameters respectively. The electric field value is calculated at the center of the cell edge, and the equivalent conductivity is obtained by the weighted average calculation of the length occupied by different media on the corresponding edge.

9. A finite difference time domain modeling system for undulating terrain and irregular anomalies, characterized in that: The finite difference time domain modeling method for undulating terrain and irregular anomaly bodies according to any one of claims 1 to 8 comprises: The triangular face element generation unit is configured to: generate a surface model according to the point cloud data of the undulating terrain and the irregular anomaly, and divide the surface model using a triangular mesh to obtain a plurality of triangular face elements; The position information generating unit is configured to: modify the normal direction of each triangular face element to obtain the triangular face element position information of the undulating terrain and the irregular abnormal body respectively; The normal vector calculation unit is configured to: import the position information of the triangular surface element into the FDTD three-dimensional forward model, and calculate the normal vector of each triangular surface element; The ray tracing unit is configured to: use FDTD grid to divide the calculation model, the FDTD grid is regarded as a ray shooting from one side of the model to the other side along the three coordinate axes of x, y and z, record the starting point and direction vector of each ray, and use the minimum storage ray and triangle intersection algorithm to obtain the intersection point of each triangle face element and the ray; The entry and exit attribute judgment unit is configured to: determine the entry and exit attributes of all intersection points according to the inner product of the normal vector of each triangular face element and the ray; The equivalent conductivity calculation unit is configured to: calculate the equivalent conductivity at the edges of the conformal grid using the conformal grid according to the in and out attributes of all intersection points; The modeling result generating unit is configured to: substitute the equivalent conductivity into the electric field iteration equation of the FDTD algorithm, perform transient electromagnetic three-dimensional forward calculation, and obtain the FDTD modeling result.

Images