Stable solution method of magnetotelluric field based on exponential scaling and renormalization

Abstract

The application relates to a magnetotelluric field stable solving method based on exponential scaling and renormalization, and belongs to the technical field of geophysics. The method comprises the following steps: obtaining three-dimensional conductivity tensors of each layer and constructing two-dimensional equivalent horizontal conductivity tensors; determining an exponential growth characteristic quantity based on layer thickness and the two-dimensional tensors, and scaling an original field propagation matrix to generate stable layer matrices of each layer according to the exponential growth characteristic quantity; multiplying the stable layer matrices and lower layer propagation matrices layer by layer from the bottom layer of a model upwards, normalizing an obtained intermediate matrix, and obtaining normalized propagation matrices of each layer interface and cumulative logarithmic scales; determining decay mode coefficients in a basement half space by using a surface electric field boundary condition and a basement decay condition, combining the normalized propagation matrix and the logarithmic scale of the topmost interface, and finally reconstructing horizontal electromagnetic field components at each layer interface based on the coefficients, the normalized propagation matrices of each interface and the logarithmic scales. The method can stably solve electromagnetic fields, scale and normalize layer by layer, and faithfully output deep electric fields.

Description

Technical Field

[0001] This application relates to the technical field of geophysics, and in particular to a method for solving the stability of the magnetotelluric field based on exponential scaling and renormalization. Background Technology

[0002] Magnetotellurics (MT) is an important geophysical exploration method that uses observations of natural alternating electromagnetic fields at the Earth's surface to detect the electrical structure of subsurface media. With the continuous improvement of exploration accuracy, three-dimensional MT forward modeling has become an important means of interpreting measured data.

[0003] Three-dimensional magnetotelluric forward modeling requires the application of boundary conditions that satisfy physical laws at the boundaries of the computational domain. Typically, the medium near the lateral boundaries is approximated as a one-dimensional layered model. After giving the polarization mode of the surface plane wave, the horizontal electric field at all depths along each boundary pillar is calculated and written into the boundary edges of the three-dimensional staggered grid. These boundary electric fields not only directly limit the outer boundary values, but also form an equivalent excitation source that drives the solution of the internal field through interaction with the three-dimensional discrete operator. Therefore, the calculation accuracy and stability of the one-dimensional horizontal electromagnetic field profile directly affect the overall reliability of the three-dimensional forward modeling.

[0004] In a one-dimensional anisotropic layered medium, the horizontal electric field and the horizontal magnetic field couple to form two independent propagation modes. The field propagation relationship of each layer can be described by a 4×4 layer transfer matrix containing hyperbolic cosine and hyperbolic sine functions. Pek and Santos (2002) established a classical method for field component propagation, which can simultaneously obtain complete information on the horizontal electric field and the horizontal magnetic field. However, this method faces a serious numerical overflow problem. On the one hand, the hyperbolic function in the layer transfer matrix contains positive exponential terms. When the frequency is high, the layer thickness is large, or the conductivity is high, the value of the positive exponential term far exceeds the range of double-precision floating-point representation. Directly calculating the hyperbolic function will produce infinity or non-numerical values, causing the single-layer matrix to become invalid before it even participates in the recursion. On the other hand, even if each single-layer matrix is ​​within the representable range, the magnitude of the elements may increase exponentially with the number of layers when the multiplication of multiple layers is continuous, again exceeding the representable range, which is particularly prominent for practical multilayer models.

[0005] To address the aforementioned issues, Pek and Santos, in their impedance recursion method, stabilized the calculation of surface impedance by simultaneously eliminating common positive exponential factors from both the numerator and denominator. However, the output of this strategy is only the surface impedance tensor, failing to provide the explicit electric field profiles at various depths required for 3D MT boundary conditions. If their field matrix propagation scheme is directly adopted, the two-stage overflow problem persists. Therefore, a method is urgently needed that can maintain the original field propagation relationship while resolving both single-layer matrix construction overflow and multi-layer matrix accumulation overflow problems, and stably output the horizontal electromagnetic field components at various depths to directly serve the boundary condition generation for 3D MT forward modeling. Summary of the Invention

[0006] This application addresses, to at least some extent, one of the technical problems in the related art.

[0007] Therefore, this application aims to provide a magnetotelluric field stability solution method based on exponential scaling and renormalization. By removing common positive exponential factors before the formation of the layer matrix to eliminate single-layer hyperbolic function overflow, implementing renormalization and recording logarithmic scale during the layer-by-layer accumulation process to suppress multi-layer product overflow, and using relative logarithmic scale to recover the horizontal electromagnetic field components at each depth, the method achieves numerical stability calculation of the electromagnetic field transfer along the entire path from the basement to the surface, and directly outputs the boundary electric field profile for three-dimensional magnetotelluric forward modeling.

[0008] To achieve the above objectives, this application provides a method for solving the stability of the magnetotelluric field based on exponential scaling and renormalization, including: Obtain the three-dimensional conductivity tensor and thickness of each layer in the layered medium model, and construct the two-dimensional equivalent horizontal conductivity tensor of each layer based on the three-dimensional conductivity tensor of each layer. Based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, the exponential growth characteristic of the corresponding layer is determined, and the original field propagation matrix of the corresponding layer is scaled based on the exponential growth characteristic to generate the stabilization layer matrix of each layer. Starting from the bottom layer of the layered medium model, the intermediate matrix corresponding to the current layer is obtained by going up layer by layer based on the stabilization layer matrix of the current layer and the propagation matrix accumulated by the lower layer. The intermediate matrix corresponding to the current layer is normalized and the corresponding logarithmic scale information is recorded to obtain the normalized propagation matrix of each layer interface and its corresponding cumulative logarithmic scale. Based on the given surface electric field boundary conditions and the field attenuation conditions of the base half-space, the attenuation mode coefficients in the base half-space are obtained using the normalized propagation matrix of the topmost interface and the logarithmic scale information. Based on the attenuation mode coefficients, the normalized propagation matrices of each interface, and the logarithmic scaling information, the horizontal electromagnetic field components at each interface are obtained.

[0009] In the technical solution, a two-dimensional equivalent horizontal conductivity tensor is first constructed using the three-dimensional conductivity tensors of each layer. This eliminates redundant vertical electrical parameters unrelated to one-dimensional electromagnetic field propagation, simplifying subsequent modal and propagation constant calculations while fully preserving the anisotropic response characteristics of the formation. Then, based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, exponentially growing characteristic quantities are determined, and the original field propagation matrix is ​​scaled to generate a stable layer matrix. This allows for early correction of exponential terms during the single-layer matrix construction stage, preventing floating-point overflow caused by direct calculation of hyperbolic functions, and adjusting only the overall matrix amplitude without altering the coupling relationship of field components. Subsequently, the stable layer matrix is ​​multiplied by the lower-layer propagation matrix from the bottom up to obtain an intermediate matrix, which is simultaneously normalized and its logarithmic scale is recorded. This continuously suppresses the exponential amplification of amplitude caused by continuous multiplication of multiple matrices, retaining all scale information uniformly using scalar logarithms and reducing memory usage. Meanwhile, the relative proportions of the field components within the matrix remain unchanged. Next, combining the surface electric field boundary and the basement attenuation condition, the basement attenuation mode coefficients are solved using the top-level normalized propagation matrix and logarithmic scale. The surface electric field required for the three-dimensional simulation is used as input, and the global scale is integrated into the coefficients to be solved, ensuring that the numerical range of the matrix involved in the calculation is smooth and reducing the ill-conditionedness that occurs when inverting the linear equation system matrix. Finally, the horizontal electromagnetic field components of each interface are recovered by combining the mode coefficients, the normalized propagation matrix of each layer, and the logarithmic scale. The ultra-large exponential calculation is avoided by relying on the relative logarithmic scale between the interfaces. The electromagnetic fields of each depth output from top to bottom automatically satisfy the boundary conditions at the interlayer interfaces. The obtained horizontal electric field can be directly used for three-dimensional mesh boundary assignment without interpolation transformation. While achieving numerical stability throughout the process, the depth electric field profile required for three-dimensional magnetotelluric simulation is completely output.

[0010] In some embodiments of this application, determining the exponentially growing characteristic quantities of different electromagnetic modes in the corresponding layer based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor includes: The two-dimensional equivalent horizontal conductivity tensor of the corresponding layer is subjected to eigenvalue decomposition to obtain the principal conductivity eigenvalues ​​corresponding to the two modes of the corresponding layer, wherein the two modes include the first mode and the second mode. Based on the two principal conductivity eigenvalues, determine the complex vertical propagation constants corresponding to the two principal conductivity eigenvalues ​​of the corresponding layer; Multiply the two complex vertical propagation constants by the thickness of the corresponding layer and take the real part to obtain the corresponding exponential growth characteristic.

[0011] In the technical solution, two independent propagation modes are separated by intrinsic decomposition of the two-dimensional equivalent horizontal conductivity tensor. The corresponding complex vertical propagation constant is accurately solved, and exponential growth-related feature values ​​are extracted in combination with the layer thickness. After offsetting, the calculated values ​​of all modal exponents can fall within the floating-point safe range, which strictly conforms to the real law of dual-mode electromagnetic propagation in the formation. The entire layer's numerical preprocessing can be completed by storing only a single scalar offset parameter. While achieving pre-suppression of single-layer exponential overflow, the program's memory read and write overhead is reduced, and the anisotropic electromagnetic response characteristics of the medium itself are not destroyed.

[0012] In some embodiments of this application, scaling the original field propagation matrix of the corresponding layer based on the exponentially growing feature quantity to generate a stable layer matrix for each layer includes: Determine the original hyperbolic cosine function value and the original hyperbolic sine function value corresponding to each mode respectively; Using the exponentially growing feature, the original hyperbolic cosine function value and the original hyperbolic sine function value of each mode are exponentially scaled to obtain the scaled hyperbolic cosine function value and the scaled hyperbolic sine function value. Based on the complex vertical propagation constant and the corresponding principal conductivity eigenvalue of each mode, the characteristic impedance of each mode is determined respectively; Based on the scaled hyperbolic cosine function value, the scaled hyperbolic sine function value, and the characteristic impedance of each mode, a scaling layer matrix of the corresponding layer in the local principal axis coordinate system is constructed. Transform the scaling layer matrix to the global coordinate system to obtain the stabilization layer matrix of the corresponding layer.

[0013] In the technical solution, the parameters of the original hyperbolic function are first calculated by combining the complex vertical propagation constants of each mode with the layer thickness. This accurately matches the basic mathematical form of the analytical solution for electromagnetic propagation in the strata. By relying on the exponentially growing characteristic quantity, the corrected hyperbolic function is obtained through exponential scaling. This directly avoids the floating-point over-limit problem caused by the excessively large original exponent. The characteristic impedance is solved by combining the principal conductivity and the propagation constant, which can completely restore the inherent coupling ratio between the electric field and the magnetic field. The scaled layer matrix in the local principal axis coordinate system is constructed using the numerically corrected hyperbolic function and the characteristic impedance. This not only fully preserves the physical interaction relationship of the cross-coupling of the two sets of modes, but also controls the numerical range of the matrix elements throughout the process. Finally, the matrix in the principal axis is mapped to the global coordinate system through coordinate transformation, resulting in a stabilized layer matrix that can be directly adapted to the calculation of any surface polarized electric field. Only the numerical amplitude is uniformly corrected throughout the process, without changing the propagation law of the electromagnetic field itself.

[0014] In some embodiments of this application, constructing the scaling layer matrix of the corresponding layer in the local principal axis coordinate system based on the scaled hyperbolic cosine function value, the scaled hyperbolic sine function value, and the characteristic impedance of each mode includes: Based on the scaled hyperbolic cosine function value corresponding to the first mode, the scaled hyperbolic sine function value corresponding to the first mode, the scaled hyperbolic cosine function value corresponding to the second mode, and the scaled hyperbolic sine function value corresponding to the second mode, as well as the characteristic impedance and the inverse of the characteristic impedance, two sets of two-dimensional field transfer sub-matrices are constructed. The two sets of two-dimensional field transfer sub-matrices are merged to obtain the scaling layer matrix of the corresponding layer in the local principal axis coordinate system.

[0015] In the technical solution, two sets of two-dimensional field transfer sub-matrices are constructed by combining the scaled hyperbolic cosine and hyperbolic sine function values, as well as the characteristic impedance and its inverse. The sub-matrices are then merged to obtain the scaled layer matrix. This not only fully adapts to the field propagation characteristics of the two sets of modes decoupled in the local principal axis coordinate system, but also avoids the numerical overflow problem caused by directly solving the exponential field general solution, ensuring the numerical calculation stability of the layer matrix. At the same time, it realizes the construction of a 4th-order scaled layer matrix. The logic is clear and easy to implement in the program. It can also accurately characterize the field quantity transfer relationship of electromagnetic waves from the bottom to the top of the dielectric layer.

[0016] In some embodiments of this application, the step of exponentially scaling the original hyperbolic cosine function values ​​and the original hyperbolic sine function values ​​of each mode using the exponentially growing feature to obtain scaled hyperbolic cosine function values ​​and scaled hyperbolic sine function values ​​includes: Compare the magnitudes of the two exponentially growing feature quantities in the corresponding layer, and select the maximum value as the common offset of the corresponding layer; For each mode in each layer, the complex vertical propagation constant of that mode is multiplied by the thickness of that layer to obtain the complex propagation thickness value of that mode; Construct positive exponent terms with the complex propagation thickness value as the exponent, and negative exponent terms with the negative value of the complex propagation thickness value as the exponent, respectively; Subtract the common offset from the exponent portion of the positive exponent term and the exponent portion of the negative exponent term respectively to obtain the scaled positive exponent term and the scaled negative exponent term. The hyperbolic cosine function value of this mode is obtained based on the sum of the scaled positive exponent term and the scaled negative exponent term. The modally scaled hyperbolic sine function value is obtained based on the difference between the scaled positive exponent term and the scaled negative exponent term.

[0017] In the technical solution, the maximum values ​​of two sets of eigenvalues ​​are selected as a unified exponential offset benchmark. This can completely preserve the imaginary part of the electromagnetic field phase correlation and only constrain the growing real part. The exponential growth eigenvalue is uniformly set as the common offset of this layer, which can achieve the same set of correction parameters to adapt to the two propagation modes within the layer. There is no need to set correction parameters separately, which reduces storage overhead. The product of the propagation constant of each mode and the layer thickness is calculated first to obtain the complex propagation thickness, which can accurately describe the complete propagation process of the electromagnetic field in the stratum. Positive and negative exponential terms are constructed to completely restore the original mathematical structure of the hyperbolic function. The common offset is directly subtracted from the exponential part instead of the result of the scaling function after the fact. This avoids the risk of numerical failure caused by the original super-large exponent being calculated first. After the offset, the real parts of the positive and negative exponential terms are not greater than 0 and fall within the stable calculation range of double-precision floating point. Then, the scaled hyperbolic cosine and hyperbolic sine are generated by summing and subtracting the exponential terms after the offset. The resulting functions differ from the original hyperbolic function by a fixed constant multiple. The relative amplitude ratio of each component of the electromagnetic field does not change, and the true field response characteristics are completely preserved.

[0018] In some embodiments of this application, the intermediate matrix corresponding to the current layer is normalized, and the corresponding logarithmic scale information is recorded, including: Extract the maximum value among all elements in the intermediate matrix and use it as the normalization factor for this normalization process; Divide each element of the intermediate matrix by the normalization factor to obtain the normalized propagation matrix of the current layer interface; The natural logarithm of the normalization factor and the common offset corresponding to the current layer are added to the cumulative logarithmic scale accumulated in the lower layer to obtain the logarithmic scale information of the current layer interface.

[0019] In the technical solution, the maximum magnitude of all elements of the intermediate matrix is ​​extracted as the normalization factor, which can accurately lock the overall magnitude of the matrix. After dividing all elements of the matrix synchronously by this factor, the maximum magnitude of the matrix can be constrained to 1, which continuously suppresses the problem of continuous amplification of magnitude during the multiplication of multiple matrices. The logarithmic value of the normalization factor is superimposed with the common offset of the current layer and incorporated into the existing cumulative logarithmic scale. The entire layer scale change information is uniformly contained by a single numerical variable. Compared with storing multiple sets of scaling matrices, the memory usage is greatly reduced. The logarithmic scale only performs linear accumulation operation and does not have its own overflow risk. Moreover, the entire matrix is ​​scaled synchronously only, and the relative proportions of each electromagnetic field component inside the matrix will not produce distortion deviation.

[0020] In some embodiments of this application, obtaining the attenuation mode coefficients in the base half-space based on the given surface electric field boundary conditions and the field attenuation conditions of the base half-space, using the normalized propagation matrix of the topmost interface and the logarithmic scaling information, includes: Based on the electrical parameters of the substrate half-space, construct the substrate field matrix that retains only the two downward decaying electromagnetic modes in the substrate half-space; Using the normalized propagation matrix of the topmost interface, the logarithmic scaling information of the topmost interface, and the base field matrix, a mapping matrix from the attenuation mode coefficients to the horizontal electric field of the ground surface is constructed. The exponential factor of the cumulative logarithmic scale of the topmost interface is incorporated into the decay mode coefficients to be solved. The attenuation mode coefficients are obtained by solving the linear equations formed by the mapping matrix based on the given surface electric field boundary conditions.

[0021] In the technical solution, a base field matrix is ​​constructed based on the base electrical parameters, retaining only the downward decaying modes. This conforms to the physical law of natural decay of electromagnetic fields in deep strata, avoiding interference from growth modes that have no physical meaning. By combining the top-level normalized propagation matrix, logarithmic scale, and base field matrix, a mapping relationship between modes and the surface electric field is established. The global scale exponential factor is integrated into the coefficients of the modes to be determined, so that the elements of the mapping matrix maintain a moderate numerical range and reduce the ill-conditionedness that occurs when inverting the linear equation system matrix. The base mode coefficients can be directly output by solving the equation system based on the standard surface polarization electric field conditions. The dimension of the equation system is only second order, which greatly reduces the computational cost of solving the equations.

[0022] In some embodiments of this application, obtaining the horizontal electromagnetic field components at each interface based on the attenuation mode coefficients, the normalized propagation matrix of each interface, and the logarithmic scaling information includes: For the target interface, the field state without scaling is calculated using the normalized propagation matrix corresponding to the interface, the basis field matrix, and the attenuation mode coefficients. The relative scale factor is determined based on the difference between the cumulative logarithmic scale corresponding to the interface and the cumulative logarithmic scale corresponding to the surface interface. The field state before scaling is scaled using the relative scale factor to obtain the horizontal electromagnetic field component at the interface.

[0023] In the technical solution, the original field state without scale correction can be restored by taking the normalized propagation matrix and the base field matrix corresponding to the target interface and combining them with the attenuation mode coefficient. The relative scale factor is constructed by using the difference between the two sets of cumulative logarithms of the target interface and the surface interface. The problem of exceeding the limit caused by global large-scale numerical calculation can be avoided by relying solely on the difference to participate in the exponential operation. After combining the unscaled field state with the relative scale factor for correction, the complete horizontal electromagnetic field components of the interface can be obtained directly. The single-depth field solution is completed independently by relying on the logarithmic parameters stored independently by each interface. The output electromagnetic field naturally satisfies the physical condition of interlayer continuity. The obtained horizontal electric field can be directly used for the layout of the three-dimensional mesh boundary without additional interpolation, conversion and other post-processing steps, simplifying the calculation process.

[0024] In some embodiments of this application, constructing the two-dimensional equivalent horizontal conductivity tensor of each layer based on the three-dimensional conductivity tensor of each layer includes: The vertical conductivity component, horizontal conductivity component, and coupled conductivity component between the vertical and horizontal components are determined based on the three-dimensional conductivity tensor. Based on the characteristic that the electromagnetic field in a one-dimensional layered medium remains unchanged along the horizontal direction, and using the physical constraint that the vertical current is zero, the correlation between the vertical electric field and the horizontal electric field is established. Based on the aforementioned correlation, the vertical electric field component of the three-dimensional conductivity tensor is eliminated to obtain a two-dimensional equivalent horizontal conductivity tensor containing only the horizontal component. Each component of the two-dimensional equivalent horizontal conductivity tensor is obtained by subtracting a correction term from the horizontal conductivity component. The correction term is the quotient obtained by dividing the product of the coupled conductivity components by the vertical conductivity component.

[0025] In the technical solution, the horizontal, vertical, and cross-coupled electrical components within the three-dimensional conductivity tensor are first completely decomposed to accurately decompose the omnidirectional conductivity characteristics of the formation. Combined with the inherent physical constraints of the one-dimensional model where the electromagnetic field remains unchanged horizontally and the vertical current is zero, an electric field correlation equation is established. This completes the elimination and simplification of the vertical electric field components within the three-dimensional tensor, removing vertically related electrical parameters that do not play a role in the one-dimensional propagation process. The simplified result is a two-dimensional equivalent horizontal conductivity tensor that only acts on the horizontal electric field. Each component of the tensor is calculated by subtracting the correction term composed of the coupling term and the vertical conductivity from the original horizontal conductivity. This fully preserves the anisotropic cross-conductivity characteristics of the formation, while compressing the computational scale of subsequent mode decomposition and propagation constant solutions. It also ensures that the simplified tensor maintains the mathematical properties of symmetry and positive definiteness, avoiding abnormal results such as non-physical negative conductivity and complex eigenvalues ​​in subsequent solutions.

[0026] In some embodiments of this application, the solution method further includes: The horizontal electromagnetic field components at each interface are written into the top boundary edge, bottom boundary edge, and corresponding sidewall boundary edge of the three-dimensional staggered grid according to the depth position to form a known boundary field vector. Based on the known boundary field vector, the preset three-dimensional internal field equation is solved to obtain the internal electric field under each polarization direction; The known boundary field vector is superimposed with the internal electric field to obtain the electric field distribution inside the three-dimensional computational domain.

[0027] In the technical solution, the corresponding horizontal electromagnetic fields at each depth are filled into the various boundary edges of the three-dimensional staggered grid and the boundary field vectors are assembled. This can directly match the boundary constraint format of the three-dimensional magnetotelluric discrete grid without additional conversion or interpolation processing. Based on this boundary field, the excitation terms of the three-dimensional equation are constructed to solve the unknown internal electric field. This can achieve physical self-consistency between the boundary constraints and the global field control equations. Finally, by superimposing the boundary field and the solved internal electric field, the complete global electric field distribution can be obtained. This fully realizes the entire process from one-dimensional stable field calculation to three-dimensional forward simulation. The output electric field data can be directly used for the extraction of electrical parameters such as apparent resistivity and impedance, and for the interpretation of geological structures.

[0028] Compared to related technologies, the magnetotelluric field stabilization solution method based on exponential scaling and renormalization provided in this application first constructs a two-dimensional equivalent horizontal conductivity tensor using the three-dimensional conductivity tensor of each layer. This eliminates redundant vertical electrical parameters unrelated to one-dimensional electromagnetic field propagation, simplifying subsequent mode and propagation constant calculations while fully preserving the anisotropic response characteristics of the strata. Then, based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, it determines the exponential growth characteristic quantity and scales the original field propagation matrix to generate a stabilization layer matrix. This allows for early correction of the exponential term during the single-layer matrix construction stage, preventing floating-point overflow caused by direct calculation of hyperbolic functions, and only adjusting the overall matrix amplitude without changing the coupling relationship of field components. Subsequently, the stabilization layer matrix is ​​multiplied by the lower-layer propagation matrix from the bottom up to obtain an intermediate matrix, which is then simultaneously normalized and its logarithmic scale recorded. This continuously suppresses the exponential amplification of amplitude caused by continuous multiplication of multiple matrices, achieving scalar-logarithmic unification. The process retains all scale information, reduces memory usage, and maintains the relative proportions of field components within the matrix. Then, combining the surface electric field boundary and basement attenuation conditions, the basement attenuation mode coefficients are solved using the top-level normalized propagation matrix and logarithmic scale. The surface electric field required for 3D simulation is used as input, and the global scale is integrated into the coefficients to be calculated, ensuring a smooth range of matrix values ​​and reducing ill-conditioned phenomena during matrix inversion of linear equations. Finally, the horizontal electromagnetic field components at each interface are recovered by combining the mode coefficients, the normalized propagation matrices of each layer, and the logarithmic scale. The large exponential calculations are avoided by relying on the relative logarithmic scale between interfaces. The electromagnetic fields at each depth output from top to bottom automatically satisfy the boundary conditions at the interlayer interfaces. The resulting horizontal electric field can be directly used for 3D mesh boundary assignment without interpolation transformation, achieving full numerical stability throughout the process while completely outputting the depth electric field profile required for 3D magnetotelluric simulation.

[0029] As can be seen from the above technical solutions, additional aspects and advantages of this application will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of this application. Attached Figure Description

[0030] Figure 1 This is a schematic diagram of a one-dimensional layered medium model with interface depth, finite layer thickness and equivalent horizontal conductivity tensor of each layer according to an embodiment of this application. Figure 2 This is a flowchart illustrating the method for solving the stability of the magnetotelluric field according to the embodiments of this application; Figure 3 In this embodiment, the COMMEMI-3D-2 model is used. Under the condition of a frequency of 0.0001Hz, along the x=0 survey line (the horizontal axis is the y-coordinate), the method of this application is combined with the EM3DANI method to obtain... A schematic diagram showing the results of comparing the resistivity of the modes; Figure 4The embodiments of this application employ COMMEMI 3D 2. Model: Under the condition of a frequency of 0.0001Hz, along the survey line x=0 (the horizontal axis is the y-coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the phase of the mode impedance; Figure 5 In this embodiment, the COMMEMI-3D-2 model is used. Under the condition of a frequency of 0.01Hz, along the x=0 survey line (the horizontal axis is the y-coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the resistivity of the modes; Figure 6 The embodiments of this application employ COMMEMI 3D 2. The model, under the condition of a frequency of 0.01Hz, along the survey line x=0 (the horizontal axis is the y-coordinate), combines the method of this application with the EM3DANI method. A schematic diagram showing the results of comparing the phase of the mode impedance; Figure 7 In this embodiment, the COMMEMI-3D-2 model is used. Under the condition of a frequency of 1Hz, along the x=0 measurement line (the horizontal axis is the y-coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the resistivity of the modes; Figure 8 The embodiments of this application employ COMMEMI 3D 2. Model: Under the condition of 1Hz frequency, along the x=0 survey line (the horizontal axis is the y-coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the phase of the mode impedance; Figure 9 In this embodiment, the COMMEMI-3D-2 model is used. Under the condition of a frequency of 100Hz, along the x=0 survey line (the horizontal axis is the y-coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the resistivity of the modes; Figure 10 The embodiments of this application employ COMMEMI 3D 2. Model: Under the condition of 100Hz, along the survey line x=0 (horizontal axis is y coordinate), the method of this application and the EM3DANI method are compared. A schematic diagram showing the results of comparing the mode impedance phase. Detailed Implementation

[0031] In the description of this application, it should be understood that the terms "center", "longitudinal", "lateral", "length", "width", "thickness", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", "clockwise", "counterclockwise", "axial", "radial", "circumferential", etc., indicating the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, are only for the convenience of describing this application and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of this application.

[0032] In this application, unless otherwise expressly specified and limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral part; they can refer to a mechanical connection, an electrical connection, or a connection that allows communication between components; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication between two components or the interaction between two components, unless otherwise expressly limited. Those skilled in the art can understand the specific meaning of the above terms in this application based on the specific circumstances.

[0033] In this application, unless otherwise expressly specified and limited, "above" or "below" the second feature can mean that the first feature is in direct contact with the second feature, or that the first feature is in indirect contact with the second feature through an intermediate medium. Furthermore, "above," "on top of," and "over" the second feature can mean that the first feature is directly above or diagonally above the second feature, or simply that the first feature is at a higher horizontal level than the second feature. "Below," "below," and "under" the second feature can mean that the first feature is directly below or diagonally below the second feature, or simply that the first feature is at a lower horizontal level than the second feature.

[0034] In this application, the terms "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., refer to a specific feature, structure, material, or characteristic described in connection with that embodiment or example, which is included in at least one embodiment or example of this application. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0035] The present application will now be described in detail through exemplary embodiments. However, it should be understood that, without further description, elements, structures, and features in one embodiment may be advantageously incorporated into other embodiments.

[0036] It should be noted that in the passive conductive region, the time factor is taken. The frequency domain magnetotelluric field control equation is:

[0037] in, For the complex electric field, For curl operator, The permeability of free space, Angular frequency, Let be the conductivity tensor. The imaginary unit, .

[0038] Discretizing this equation yields a set of three-dimensional discrete electric field operator equations:

[0039] in, For a three-dimensional discrete electric field operator, Let be a complex vector consisting of all the degrees of freedom of the electric field.

[0040] The total electric field over the entire domain is divided into two parts: a known boundary field and an unknown interior field. These two parts satisfy the following:

[0041] in It exists only within the computational domain and has a value of zero at its outer boundary. With electric field values ​​existing only at the grid boundary edges and zero in all other degrees of freedom, the internal field solution equations can be derived by substituting the decomposition formula into the discretization equations:

[0042] This clarifies that the boundary electric field not only directly defines the outer boundary field value of the 3D mesh, but also couples with the discrete operator to generate an equivalent right-hand side term that drives the solution of the unknown internal electric field. In 3D magnetotelluric numerical simulations, the sidewalls and top / bottom boundary media of the computational domain are usually approximated as a one-dimensional columnar model whose electrical properties vary only with depth. Given the polarization conditions of the surface plane waves, the entire domain depth must be solved along each boundary one-dimensional column. Only with a horizontal electric field profile can the boundary values ​​of the three-dimensional mesh be assigned. The computational stability and output completeness of the one-dimensional field profile directly determine the accuracy of the three-dimensional forward modeling.

[0043] Pek and Santos established a classical method for field component propagation, which can simultaneously obtain complete information on both the horizontal electric and magnetic fields. Existing classical Pek and Santos one-dimensional anisotropic magnetotelluric theory consists of two systems: impedance recursion and four-component field propagation. Impedance recursion can only stably solve for surface impedance parameters and cannot output depth-varying parameters. The electric field profile, and the four-component field transfer matrix completely preserves all depth information of the electromagnetic field, naturally adapts to the requirements of three-dimensional Dirichlet boundary assignment, and is the basic theoretical framework supporting multidimensional MT simulation. Under this framework, the field matrix of each layer is constructed by hyperbolic cosine and hyperbolic sine functions. During the propagation of the electromagnetic field along the depth of the stratum, two sets of exponential modes, growth and decay, will be generated simultaneously.

[0044] In existing technologies, numerical calculations using the original four-component field matrices of Pek and Santos suffer from two levels of floating-point overflow defects, and the real part of the propagation constant is also affected under high-frequency, thick, and high-conductivity strata conditions. The value is too large. Breaking the upper limit of IEEE double-precision floating-point expression, the single-layer matrix construction stage generates infinite and non-numerical illegal results; even if the single-layer matrix value is within the legal range, the amplitude of matrix elements will expand exponentially during the continuous multiplication of multiple 4×4 matrices, resulting in numerical overflow again; Pek and Santos' exponential elimination stabilization method for impedance fractions is only applicable to the synchronous simplification of impedance numerator and denominator, and cannot be transferred to the entire process of constructing explicit four-component field matrices, layer-by-layer accumulation, and complete recovery of the depth field. Moreover, the existing stabilization algorithm can only output the surface impedance and cannot generate the entire depth horizontal electric field required for the three-dimensional mesh, and cannot completely realize the technical link from one-dimensional field calculation to three-dimensional magnetotelluric forward modeling simulation.

[0045] Based on this, this application proposes a magnetotelluric field stability solution method based on exponential scaling and renormalization. By eliminating the vertically correlated electrical components of the three-dimensional symmetric positive definite conductivity tensor to obtain a two-dimensional equivalent horizontal tensor and completing mode decomposition, the maximum real part exponents of two propagation modes are extracted as common offsets before each layer matrix is ​​constructed. The scaling hyperbolic function is directly constructed to form a stable layer matrix. After performing layer matrix accumulation from the base upwards, real-time renormalization is performed using the maximum modulus of the matrix elements and the logarithmic scale is uniformly stored. The base attenuation mode is solved by combining the given surface electric field. The full-depth four-component electromagnetic field is restored by relying on the relative logarithmic scale between the interfaces. Finally, the depth electric field is written into the three-dimensional staggered grid to generate boundary vectors and equivalent right-hand terms. This achieves simultaneous suppression of numerical overflow at both single-layer and multi-layer matrix levels while completely preserving the original electromagnetic field physical transmission relationship. This solves the technical problems of traditional field matrix calculations exceeding the floating-point limit at both levels and only being able to output the surface impedance, which cannot supply the full-domain depth boundary electric field of the three-dimensional magnetotelluric field.

[0046] In the following, embodiments of this application will be described in detail with reference to the accompanying drawings.

[0047] As attached Figure 2 As shown in the illustrative embodiment of a magnetotelluric field stability solution method based on exponential scaling and renormalization in this application, the solution method includes: Step S100: Obtain the three-dimensional conductivity tensor and thickness of each layer in the layered medium model, and construct the two-dimensional equivalent horizontal conductivity tensor of each layer based on the three-dimensional conductivity tensor of each layer.

[0048] It should be noted that, referring to Figure 1 A layered medium model, also known as a one-dimensional layered model, refers to a stratigraphic model in which the electrical parameters of the subsurface medium change only along the depth direction, while the electrical properties remain uniform in the horizontal direction. To avoid confusion regarding layer numbers, interface numbers, mode numbers, and propagation directions, this invention adopts a unified notation convention, defining the Earth's surface as... The z-axis points vertically downwards, as shown below. Figure 1 As shown, this one-dimensional layered model includes A layer of finite thickness and a uniform substrate half-space beneath it.

[0049] For the first The thickness of a finite layer, For the first The depth of the interface, where the first... The numbering of each finite layer is taken The top interface of this layer is The bottom interface is The thickness meets the requirements. The two-dimensional equivalent horizontal conductivity tensor of this layer is .

[0050] The base half-space is numbered as Its top position is The interface numbering adopts a unified format. ,in Represents the Earth's surface.

[0051] Specifically, the vertical conductivity component, horizontal conductivity component, and vertical-horizontal coupled conductivity component are first extracted based on the three-dimensional conductivity tensor. Then, based on the physical constraints that the horizontal gradient of the electromagnetic field in the one-dimensional layered medium is zero and the vertical current of the formation is zero, an equation relating the vertical and horizontal electric fields is established. Next, a two-dimensional equivalent horizontal conductivity tensor is generated by eliminating vertical electric field-related parameters within the three-dimensional tensor through elimination. This step can eliminate redundant vertical electrical parameters unrelated to the propagation of the one-dimensional electromagnetic field, simplifying subsequent modal and propagation constant calculations while fully preserving the anisotropic response characteristics of the formation. The simplified two-dimensional tensor maintains symmetric positive definite properties, avoiding non-physical negative conductivity and complex eigenvalue anomalies in subsequent solutions.

[0052] In some embodiments, to achieve an equivalent conversion from a three-dimensional conductivity tensor to a two-dimensional horizontal tensor calculated using a one-dimensional layered model, and to eliminate invalid vertical electrical components during the propagation of the one-dimensional electromagnetic field, step S100, based on the three-dimensional conductivity tensors of each layer, constructs the two-dimensional equivalent horizontal conductivity tensors of each layer, including: The vertical conductivity component, horizontal conductivity component, and coupled conductivity component between the vertical and horizontal directions are determined based on the three-dimensional conductivity tensor. Based on the characteristic that the electromagnetic field in a one-dimensional layered medium remains constant along the horizontal direction, and utilizing the physical constraint that the vertical current is zero, the correlation between the vertical and horizontal electric fields is established. According to this correlation, the vertical electric field component of the three-dimensional conductivity tensor is eliminated to obtain a two-dimensional equivalent horizontal conductivity tensor containing only the horizontal component. Each component of the two-dimensional equivalent horizontal conductivity tensor is obtained by subtracting a correction term from the horizontal conductivity component; the correction term is the quotient obtained by dividing the product of the coupled conductivity components by the vertical conductivity component.

[0053] For example, for the first First, read the 3×3 symmetric positive definite three-dimensional conductivity tensor of that layer. The expression for the three-dimensional conductivity tensor is:

[0054] Extract the vertical conductivity component from the aforementioned three-dimensional conductivity tensor. Horizontal conductivity component and the coupled conductivity components between the vertical and horizontal directions Because the electromagnetic field in a one-dimensional layered model satisfies the property that the gradient along the horizontal x and y directions is zero, i.e. Electromagnetic field constraint conditions, taking time factor By combining the vertical current constraint relationship determined by Maxwell's equations in the frequency domain, a vertical electric field is established. With horizontal electric field component The relationship between them, where the frequency domain Maxwell equations simplify to:

[0055]

[0056] Among them, the positive or negative sign generated by the frequency factor is determined by the time factor. The time was agreed upon. and These represent the elimination of the vertical electric field. And the equivalent horizontal conductivity acting in the x and y directions after vertical current constraint; This represents the equivalent coupling of the electric field in the y-direction to the current in the x-direction. This represents the equivalent coupling of the electric field in the x-direction to the current in the y-direction.

[0057] Then, based on this correlation, a vertical electric field is applied to the three-dimensional conductivity tensor. Elimination process, removing all Related terms, thus yielding a 2×2 two-dimensional equivalent horizontal conductivity tensor. The calculation expressions for each component of the two-dimensional equivalent horizontal conductivity tensor are as follows:

[0058]

[0059]

[0060]

[0061] Furthermore, the horizontal electric field vector is denoted as... The horizontal magnetic field vector is denoted as Therefore, the horizontal electric field and horizontal magnetic field are combined to form a four-component horizontal electromagnetic field vector. The horizontal electromagnetic field vector satisfies After eliminating the vertical electric field, eliminating the horizontal magnetic field-related terms in the equations allows us to derive the equations relating only to the horizontal electric field. The second-order ordinary differential governing equation, which is related to the two-dimensional equivalent horizontal conductivity tensor. The matrix representation is as follows:

[0062] The above governing equations establish the constraint relationship between the variation of the horizontal electric field along the depth direction inside the layered anisotropic medium, and the two-dimensional equivalent horizontal conductivity tensor. To fully characterize the equivalent conductivity of the formation in the horizontal direction, we first completely decompose the horizontal, vertical, and cross-coupled electrical components within the three-dimensional conductivity tensor, accurately decompose the omnidirectional conductivity characteristics of the formation, and establish electric field correlation equations based on the inherent physical constraints of the one-dimensional model's electromagnetic field remaining unchanged horizontally and vertical current returning to zero. We then eliminate the vertical electrical coupling components in the three-dimensional conductivity tensor that do not affect the propagation of the electromagnetic field in the one-dimensional layered model. While preserving the anisotropic cross-conductivity characteristics of the formation, we reduce the three-dimensional electrical parameters to a two-dimensional matrix form, providing the basic matrix object for subsequent intrinsic decomposition, mode separation, and complex vertical propagation constant solution.

[0063] Furthermore, in order to filter out abnormal input parameters that do not meet the physical solution prerequisites, and to avoid problems such as calculation divergence and solution failure caused by non-positive definite matrices in subsequent matrix operations, and to ensure that the electromagnetic field solution process can proceed stably, the thickness of each layer can be verified first. Yes, all are satisfied. Based on the three-dimensional conductivity tensors of each layer, the conversion is completed, and the two-dimensional equivalent horizontal conductivity tensor is output. Then, further verification of the... Whether the positive definite condition is met is determined by the following expression: ,and If the verification fails, it indicates an anomaly in the input parameters of the current layer's three-dimensional conductivity tensor. The calculation process is then terminated, and an error message is output. In this way, the transformed two-dimensional horizontal tensor maintains symmetric positive definite properties, avoiding calculation defects such as non-physical negative conductivity and complex eigenvalues. Simultaneously, it fully preserves the cross-conductivity response characteristics caused by formation anisotropy. The pre-verification step can also identify unreasonable input parameters in advance, preventing non-physical interpretations in subsequent matrix operations and ensuring the robustness of the entire algorithm's calculation process.

[0064] Step S200: Based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, determine the exponential growth characteristic quantity of the corresponding layer, and scale the original field propagation matrix of the corresponding layer based on the exponential growth characteristic quantity to generate the stabilized layer matrix of each layer.

[0065] Based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, an exponentially growing characteristic quantity is determined, and the original field propagation matrix is ​​scaled to generate a stable layer matrix. By introducing the exponentially growing characteristic quantity to perform matrix scaling, the exponential term can be corrected in advance during the single-layer matrix construction stage, the numerical overflow phenomenon caused by the rapid growth of the exponential function can be suppressed, and the numerical distortion problem that is easy to occur when directly solving the original field propagation matrix can be avoided. The scaling operation only changes the magnitude of the matrix, and the inherent physical coupling correspondence between the electromagnetic field components in the layer is completely preserved. It does not destroy the original field propagation physical law of the layered medium, ensures the numerical stability of the subsequent layer-by-layer matrix accumulation operation, avoids the interruption of the solution due to numerical anomalies, and improves the computational reliability and result accuracy of the forward modeling solution of layered anisotropic magnetotelluric fields.

[0066] In some embodiments, determining the exponential growth characteristic quantities of different electromagnetic modes in the corresponding layer based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor in step S200 includes: The two-dimensional equivalent horizontal conductivity tensor of the corresponding layer is decomposed into eigenvalues ​​to obtain the principal conductivity eigenvalues ​​corresponding to the two modes of the corresponding layer. The two modes include the first mode and the second mode. Based on the two principal conductivity eigenvalues, the complex vertical propagation constants corresponding to the two principal conductivity eigenvalues ​​of the corresponding layer are determined. The two complex vertical propagation constants are multiplied by the thickness of the corresponding layer and the real part is taken to obtain the corresponding exponential growth characteristic quantities.

[0067] Specifically, when the first Three-dimensional conductivity tensor of the layer When it is a symmetric tensor, it satisfies , as well as Therefore, it can be deduced that the two-dimensional equivalent horizontal conductivity tensor satisfies Therefore, when the three-dimensional conductivity tensor When further defined as a symmetric positive definite tensor, the corresponding two-dimensional equivalent horizontal conductivity tensor It is also a symmetric positive definite matrix, therefore it can be used for the two-dimensional equivalent horizontal conductivity tensor. Eigenvalue decomposition is performed on the two-dimensional equivalent horizontal conductivity tensor obtained from the solution of this layer to obtain the principal conductivity eigenvalues ​​corresponding to the first and second modes of the corresponding layer. , The decomposition form is:

[0068] in, , These are two principal conductivity eigenvalues ​​that are greater than zero; Let be a unitary matrix consisting of two sets of eigenvectors, which satisfies And determinant , Physical meaning represents orbit Local principal axis rotation transformation carried out on the axis.

[0069] It should be noted that the modal numbering adopts a unified approach. , , and It is only used to characterize two main conductivity modes within the same layer, has no meaning of propagation direction, and does not represent positive or negative propagation direction.

[0070] In the You formation Within the defined principal axis coordinate system, a coordinate rotation transformation is performed on the horizontal electric field and the horizontal magnetic field to obtain:

[0071] Due to the Inner layer of the array Not with depth The change allows for decoupling, resulting in two sets of uncoupled scalar differential equations:

[0072] In the formula, the complex vertical propagation constant for:

[0073] When taking the square root of the complex vertical propagation constant, two mathematical solutions are obtained, i.e., two branches. Therefore, the real part is chosen. The mathematical branch is selected according to the physical law of the downward decay of electromagnetic fields in deep strata, avoiding the introduction of incorrect calculation results by selecting non-physical growth branches.

[0074] When the local depth within the layer is... Then the top floor Bottom layer The general solution to the two sets of uncoupled scalar differential equations is: ,

[0075] in, and The first Layer mode The complex amplitudes of the decay and growth exponent terms. Both terms must be retained in finite-thickness layers; only the basal half-space needs to be discarded as the exponent increases with depth. .

[0076] Then, by multiplying the two complex vertical propagation constants by the thickness of the corresponding layer and taking their real parts, the exponential growth characteristic quantities corresponding to each mode are obtained. and These characteristic quantities characterize the magnitude of the exponential change in field strength when electromagnetic waves pass through this layer.

[0077] In the technical solution, by performing intrinsic decomposition on the two-dimensional equivalent horizontal conductivity tensor to separate two sets of independent propagation modes, the corresponding complex vertical propagation constant is accurately solved, and the exponential growth related feature value is extracted in combination with the layer thickness. It can completely preserve the imaginary part value of the electromagnetic field phase correlation and only constrain the real part of the growth. The calculated values ​​of all modal exponents after offset can fall within the floating-point safe range, strictly conforming to the real law of dual-mode electromagnetic propagation in the formation. The entire layer numerical preprocessing can be completed by storing only a single scalar offset parameter. While achieving pre-suppression of single-layer exponential overflow, it reduces the program memory read and write overhead and does not destroy the anisotropic electromagnetic response characteristics of the medium itself.

[0078] In some embodiments, the scaling process of the original field propagation matrix of the corresponding layer based on the exponentially growing feature quantity in step S200 to generate the stabilized layer matrix of each layer includes: The original hyperbolic cosine function value and the original hyperbolic sine function value corresponding to each mode are determined respectively. Using the exponential growth eigenvalue, the original hyperbolic cosine function value and the original hyperbolic sine function value of each mode are exponentially scaled to obtain the scaled hyperbolic cosine function value and the scaled hyperbolic sine function value. Based on the complex vertical propagation constant and the corresponding principal conductivity eigenvalue of each mode, the characteristic impedance of each mode is determined respectively. Based on the scaled hyperbolic cosine function value, the scaled hyperbolic sine function value and the characteristic impedance of each mode, the scaling layer matrix of the corresponding layer in the local principal axis coordinate system is constructed. The scaling layer matrix is ​​transformed to the global coordinate system to obtain the stabilization layer matrix of the corresponding layer.

[0079] Specifically, the core calculation unit inside the original field propagation matrix is ​​the hyperbolic sine function and hyperbolic cosine function under the corresponding mode. If the hyperbolic function is directly calculated, when the layer thickness is large and the real part of the propagation constant is high, the hyperbolic function will increase exponentially, directly triggering a double-precision floating-point overflow error. Therefore, this embodiment introduces the exponentially growing characteristic quantity obtained by the above steps as a unified offset, and performs offset correction on the hyperbolic function in the exponential domain to suppress the amplitude explosion of the hyperbolic function.

[0080] First, based on the complex vertical propagation constants of the two modes of the corresponding layer. and the thickness of the layer Determine the original hyperbolic cosine function value corresponding to each mode respectively. and the original hyperbolic sine function value ,in:

[0081]

[0082] Then use the common offset The original hyperbolic cosine function values ​​and the original hyperbolic sine function values ​​for each mode are exponentially scaled to obtain scaled hyperbolic cosine function values. and scaled hyperbolic sine function values .

[0083] In some embodiments, the original hyperbolic cosine function values ​​and the original hyperbolic sine function values ​​of each mode are exponentially scaled using an exponentially growing feature, resulting in scaled hyperbolic cosine function values ​​and scaled hyperbolic sine function values, including: Compare the magnitudes of the two exponentially growing feature quantities in the corresponding layer, and select the maximum value as the common offset of the corresponding layer; for each mode in each layer, multiply the complex vertical propagation constant of the mode by the thickness of the layer to obtain the complex propagation thickness value of the mode; construct a positive exponential term with the complex propagation thickness value as the exponent, and a negative exponential term with the negative value of the complex propagation thickness value as the exponent; subtract the common offset from the exponential part of the positive exponential term and the exponential part of the negative exponential term respectively to obtain the scaled positive exponential term and the scaled negative exponential term; obtain the scaled hyperbolic cosine function value of the mode based on the sum of the scaled positive exponential term and the scaled negative exponential term; obtain the scaled hyperbolic sine function value of the mode based on the difference between the scaled positive exponential term and the scaled negative exponential term.

[0084] Specifically, the magnitudes of the two exponentially growing features in this layer are compared, and the maximum value is selected as the common offset of the corresponding layer. At the same time, ensure that the value is non-negative, i.e., define:

[0085] Then, for each mode in each layer, the complex vertical propagation constant of that mode is... With layer thickness Multiplying them yields the complex propagation thickness value for that mode. Then, positive exponential terms with the complex propagation thickness as the exponent are constructed respectively. And a negative exponential term with the negative value of the complex propagation thickness as the exponent. The common offset will be subtracted from the exponents of both the positive and negative exponent terms. We obtain the scaled positive exponent term. and the scaled negative exponent term Due to common offset The real part of the positive exponent is not less than that of the two modes. The real parts of the four scaled exponential functions are all not greater than zero, and no numerical overflow will occur.

[0086] Then, based on the sum of the scaled positive exponent terms and the scaled negative exponent terms, the scaled hyperbolic cosine function value for this mode is obtained:

[0087] Based on the difference between the scaled positive exponent term and the scaled negative exponent term, the scaled hyperbolic sine function value for this mode is obtained:

[0088] Setting the exponential growth characteristic as a common offset for this layer allows for the use of the same set of correction parameters to adapt to two propagation modes within the layer, eliminating the need for separate correction parameters and reducing storage overhead. The complex propagation thickness is obtained by first calculating the product of the propagation constant of each mode and the layer thickness, which can accurately describe the complete propagation process of the electromagnetic field within the stratum. Positive and negative exponential terms are constructed to completely restore the original mathematical structure of the hyperbolic function. The common offset is directly subtracted from the exponential part instead of the result of post-scaled function, thus avoiding the risk of numerical failure caused by calculating the original super-large exponent from the root. After offset processing, the real parts of the positive and negative exponential terms are not greater than 0, and the entire process falls within the stable calculation range of double-precision floating-point numbers. Then, the scaled hyperbolic cosine and hyperbolic sine are generated by summing and subtracting the exponential terms after offset. The resulting functions differ from the original hyperbolic function by only a fixed constant multiple, and the relative amplitude ratios of each component of the electromagnetic field do not change, thus completely preserving the true field response characteristics.

[0089] Meanwhile, based on the complex vertical propagation constant of each mode and the corresponding principal conductivity eigenvalues Determine the characteristic impedance of each mode separately:

[0090] Characteristic impedance table The proportional relationship between the electric field amplitude and the magnetic field amplitude of this mode is characterized.

[0091] In some embodiments, based on the scaled hyperbolic cosine function values, scaled hyperbolic sine function values, and characteristic impedance of each mode, a scaling layer matrix in the local principal axis coordinate system for the corresponding layer is constructed, including: Based on the scaled hyperbolic cosine function value corresponding to the first mode, the scaled hyperbolic sine function value corresponding to the first mode, the scaled hyperbolic cosine function value corresponding to the second mode, the scaled hyperbolic sine function value corresponding to the second mode, as well as the characteristic impedance and the inverse of the characteristic impedance, two sets of two-dimensional field transfer sub-matrices are constructed; the two sets of two-dimensional field transfer sub-matrices are merged to obtain the scaling layer matrix of the corresponding layer in the local principal axis coordinate system.

[0092] Specifically, based on the hyperbolic cosine function values ​​scaled for each modality. Scaled hyperbolic sine function values Characteristic impedance and the inverse of characteristic impedance Construct the scaling layer matrix of the corresponding layer in the local principal axis coordinate system.

[0093] make These are four-component principal axis field vectors. Because... After rotating to the principal axis coordinate system, the relationship between the electric and magnetic fields remains the same, that is:

[0094]

[0095] Therefore, differentiating the general solutions of the two sets of uncoupled scalar differential equations and using... ,get:

[0096]

[0097] Using the general solutions of the aforementioned electric and magnetic fields, let the top of the layer... Bottom layer Eliminate undetermined coefficients and , and utilize and From the relationship, we obtain two sets of two-dimensional field transfer sub-matrices that propagate from the bottom of the layer to the top of the layer, and their 2×2 relationship is as follows:

[0098]

[0099] These two 2×2 relationships are combined into a complete 4×4 scaling layer matrix in the local principal axis coordinate system, thus obtaining the scaling layer matrix of the corresponding layer in the local principal axis coordinate system. for:

[0100] remember Then the global field and the principal axis field satisfy the following:

[0101] Then use the coordinate transformation matrix Complete the transformation from the local principal axis coordinate system to the global coordinate system, and generate the stable layer matrix in the global coordinate system. The matrix transformation relationship satisfies:

[0102] Thus, by first combining the complex vertical propagation constants of each mode with the layer thickness to calculate the relevant parameters of the original hyperbolic function, we can accurately fit the basic mathematical form of the analytical solution for electromagnetic propagation in the strata. By relying on the exponentially growing characteristic quantity to complete the exponential scaling to obtain the corrected hyperbolic function, we can directly avoid the floating-point over-limit problem caused by the excessively large original exponential term. By combining the principal conductivity and propagation constant to solve the characteristic impedance, we can completely restore the inherent coupling ratio between the electric field and the magnetic field. Using the numerically corrected hyperbolic function and the characteristic impedance, we construct a scaled layer matrix in the local principal axis coordinate system. This not only completely preserves the physical interaction relationship of the cross-coupling of the two sets of modes, but also controls the numerical range of the matrix elements throughout the process. Finally, through coordinate transformation, we map the matrix in the principal axis to the global coordinate system to obtain a stabilized layer matrix that can be directly adapted to the calculation of any surface polarized electric field. Throughout the process, only the numerical amplitude is uniformly corrected, without changing the propagation law of the electromagnetic field itself.

[0103] Step S300: Starting from the bottom layer of the layered medium model, the intermediate matrix corresponding to the current layer is obtained by going up layer by layer based on the stabilized layer matrix of the current layer and the accumulated propagation matrix of the lower layer. The intermediate matrix corresponding to the current layer is normalized and the corresponding logarithmic scale information is recorded to obtain the normalized propagation matrix of each layer interface and its corresponding accumulated logarithmic scale.

[0104] Multiplying the stabilization layer matrix with the lower propagation matrix from the bottom up to obtain the intermediate matrix and simultaneously normalizing and recording the logarithmic scale can continuously suppress the exponential amplification of amplitude caused by continuous multiplication of multi-layer matrices. It retains all scale information in a unified scalar logarithmic manner, reduces memory usage, and maintains that the relative proportions of each field component within the matrix remain unchanged.

[0105] Specifically, definition To the base interface Propagation to the interface The physical propagation matrix, let the basis position matrix Normalized propagation matrix Initial cumulative log-scale of the base Therefore, the propagation matrix for:

[0106] For example, ,and These physical propagation matrices may overflow due to multiplication. Therefore, the program does not store them directly. Instead, it preserves the normalized propagation matrix. and real numbers , so that:

[0107] Therefore, from the bottom layer To the top The recursive calculation is performed layer by layer upwards, multiplying the current layer's stabilization matrix by the normalized propagation matrix obtained from the next layer to obtain the intermediate matrix corresponding to the current layer. This yields the complete interface. Normalized propagation matrix at the location and its corresponding cumulative logarithmic scale .

[0108] For example, for the current number The layer will stabilize the layer matrix of the current layer. Multiply by the normalized propagation matrix of the lower interface To obtain a temporary intermediate matrix :

[0109] Then, the intermediate matrix corresponding to the current layer is normalized, and the corresponding logarithmic scale information is recorded to obtain the normalized propagation matrix of each layer interface and its corresponding cumulative logarithmic scale.

[0110] In some embodiments, normalizing the intermediate matrix corresponding to the current layer and recording the corresponding logarithmic scale information includes: Extract the maximum value among all elements in the intermediate matrix as the normalization factor for this normalization process. Divide each element of the intermediate matrix by the normalization factor to obtain the normalized propagation matrix of the current layer interface. Add the natural logarithm of the normalization factor and the common offset corresponding to the current layer to the accumulated logarithmic scale of the lower layer to obtain the logarithmic scale information of the current layer interface.

[0111] Specifically, extract the intermediate matrix. The maximum value of the modulus of all 16 matrix elements is used as the normalization factor for this normalization process. The normalization factor is:

[0112] in, The matrix row and column indices, It is the maximum value among its 16 modulo values, used to construct the first... Normalization factor for matrix multiplication.

[0113] The intermediate matrix is ​​scaled using a normalization factor. Each element of the intermediate matrix is ​​then divided by the normalization factor to obtain the normalized propagation matrix of the current layer interface.

[0114] The maximum value of the modulus of all elements of the normalized propagation matrix is ​​exactly 1.

[0115] At the same time, the common offset corresponding to this layer will be... and the logarithm of the normalization factor Superimposed on the lower cumulative logarithmic scale In the middle, the cumulative logarithmic scale of the current interface is updated. :

[0116] Furthermore, the normalization factor is completed at each layer. After solving the problem, a numerical validity check is required. If a normalization factor is identified... If the result is zero or a non-finite value, the current calculation process is immediately terminated, and the corresponding layer number and frequency parameters are reported. No manually set constants are used in the calculation process to replace the matrix calculation results obtained from physical derivation.

[0117] Extracting the maximum magnitude of all elements in the intermediate matrix as a normalization factor accurately locks the overall magnitude of the matrix. Dividing all matrix elements synchronously by this factor constrains the maximum magnitude to 1, continuously suppressing the problem of amplitude amplification during multi-layer matrix multiplication. The logarithmic value of the normalization factor is superimposed with the common offset of the current layer and incorporated into the existing cumulative logarithmic scale. Relying on a single numerical variable to uniformly encompass the scale change information of the entire layer significantly reduces memory usage compared to storing multiple sets of scaling matrices. The logarithmic scale only performs linear accumulation operations, eliminating the risk of overflow. Furthermore, since only all elements of the matrix are scaled synchronously, the relative proportions between the components of the electromagnetic field remain unchanged, preventing computational distortion in the electromagnetic field response. The calculation is terminated promptly when the normalization factor becomes abnormal, preventing the abnormal matrix from being propagated downwards and causing completely erroneous electromagnetic field results in subsequent solutions.

[0118] Step S400: Based on the given surface electric field boundary conditions and the field attenuation conditions of the base half-space, the attenuation mode coefficients in the base half-space are obtained using the normalized propagation matrix and logarithmic scale information of the top interface.

[0119] After performing layer-by-layer multiplication and renormalization, the global exponential scaling factor at the topmost interface, which may be extremely large or small, will be reduced. The system is absorbed into the basis decay mode coefficients (C) to be solved, thus transforming the linear system, which might otherwise be unsolvable due to numerical overflow, into a small 2×2 linear equation system with good condition numbers, consisting only of normalized stability matrices. The solution process avoids directly calculating global exponential terms that might exceed the range of double-precision floating-point numbers. Simultaneously, the field attenuation condition of the base half-space (discarding growth terms) ensures the rationality of the underlying physical field, ultimately enabling the stable, unique, and accurate determination of the amplitude coefficients of the two attenuation modes even at extremely high frequencies or with extremely high conductivity contrast. Combining the surface electric field boundary and the base attenuation condition, the base attenuation mode coefficients are solved using the top-level normalized propagation matrix and logarithmic scaling. The surface electric field required for the 3D simulation is used as input, and the global scale is integrated into the coefficients to be determined, ensuring a smooth range of matrix values ​​involved in the calculation and reducing ill-conditioned phenomena that occur when inverting linear equations.

[0120] In some embodiments, based on given surface electric field boundary conditions and field attenuation conditions of the base half-space, the attenuation mode coefficients in the base half-space are obtained using the normalized propagation matrix and logarithmic scaling information of the topmost interface, including: Based on the electrical parameters of the base half-space, a base field matrix is ​​constructed that retains only the two downward decaying electromagnetic modes in the base half-space. Using the normalized propagation matrix of the top interface, the logarithmic scale information of the top interface, and the base field matrix, a mapping matrix from the decaying mode coefficients to the horizontal electric field of the ground surface is constructed. The exponential factor of the cumulative logarithmic scale of the top interface is absorbed into the decaying mode coefficients to be solved. The linear equations formed by the mapping matrix are solved according to the given boundary conditions of the ground surface electric field to obtain the decaying mode coefficients.

[0121] Specifically, let the local depth within the basal half-space be . These are the local coordinates of the basis half-space. Corresponding layer interface , Increased represents going deeper underground, because To ensure Limited timeframe, growth items It must be zero. Therefore, only the base half-space retains... Two downward decaying modes.

[0122] If they are The electric field amplitudes at the locations are respectively and The corresponding magnetic field components are:

[0123]

[0124] It should be noted that the positive and negative signs are determined by the relationship between the electric and magnetic fields of the two modes.

[0125] At the top of the base half-space The interface, or layer The position defines the modal amplitude vector. It is used to contain the electric field amplitudes of two downward decaying modes. and :

[0126] Substituting the amplitude vector into the modal field expression yields the complete field state vector at the interface of the lower half-space in the local coordinate system. This vector contains components of both the electric and magnetic fields, and can be written in matrix multiplication form:

[0127] Where the matrix Given the modal coupling matrix of the half-space in the local coordinate system, we can directly establish the mapping relationship between the modal amplitude and the field components in the local coordinate system:

[0128] because Defined in local depth coordinates Next, it is necessary to use a transformation matrix. Transform to the global physical coordinate system to obtain the half-space modal basis field matrix corresponding to the global coordinate system. :

[0129] Therefore, the interface in the global coordinate system The relationship between the field state vector at the top of the base, i.e., the four-component field state at the top of the base, and the attenuation mode coefficient is as follows:

[0130] Furthermore, The first column represents the unit p-mode. The corresponding field states in the principal axis coordinate system, the second column represents the unit m modes. The corresponding field state in the principal axis coordinate system; transformation matrix This is used to rotate and transform the modal field state in the principal axis coordinate system to the global physical coordinate system, obtaining the basis field matrix. This matrix retains only the two types of electromagnetic modes that decay downwards within the base half-space, discarding the growing modes that diverge with increasing depth, and satisfies the attenuation boundary condition of the field finite at infinity in the half-space.

[0131] Further define the electric field selection matrix This matrix extracts the horizontal electric field component from the complete four-dimensional field state vector, in conjunction with the overall propagation relationship of the layered medium. and ，Based on the normalized propagation matrix and logarithmic scaling information of the topmost interface, a mapping matrix from the attenuation mode coefficients to the horizontal electric field at the Earth's surface is constructed. :

[0132] Then, the accumulated logarithmic scale of the topmost interface. Exponential factors The attenuation mode coefficient C is defined as follows: (The information is absorbed into the attenuation mode coefficients to be solved.)

[0133] Finally, based on the given surface electric field boundary conditions, the system of linear equations consisting of the mapping matrix is ​​solved to obtain the attenuation mode coefficients, where the given surface electric field satisfies:

[0134] Therefore, based on Obtain the attenuation mode coefficient .when At that time, this system of equations uniquely determines the attenuation mode coefficients. .

[0135] Thus, a basefield matrix that retains only downward decaying modes is constructed based on the basefield electrical parameters, conforming to the physical law of natural decay of electromagnetic fields in deep strata, avoiding interference from growth modes without physical meaning. By combining the top-level normalized propagation matrix, logarithmic scale, and basefield matrix, a mapping relationship from modes to the surface electric field is established. The global scale exponential factor is integrated into the coefficients of the modes to be determined, so that the elements of the mapping matrix maintain a moderate numerical range and reduce the ill-conditionedness that occurs when inverting the linear equation system matrix. The basefield mode coefficients can be directly output by solving the equation system based on the standard surface polarization electric field conditions. The dimension of the equation system is only second order, which greatly reduces the computational cost of solving the equations.

[0136] Furthermore, during the equation solving process, the boundary electric field mapping matrix is ​​monitored in real time. The numerical state, if determined If the matrix exhibits a singular or near-singular ill-conditioned condition, it indicates that the stability of the equation system under the current input model and frequency parameter combination is insufficient. The calculation task for the current corresponding frequency should be terminated immediately, and the layer number, frequency, and matrix status report information should be output. Similarly, arbitrary artificial constant substitutes should not be used to understand and force a solution.

[0137] Step S500: Based on the attenuation mode coefficients, the normalized propagation matrix of each interface, and the logarithmic scale information, the horizontal electromagnetic field components at each interface are obtained.

[0138] By combining modal coefficients, normalized propagation matrices of each layer, and logarithmic scales, the horizontal electromagnetic field components of each interface are recovered. By relying on the relative logarithmic scale between interfaces to avoid ultra-large exponential calculations, the electromagnetic fields at each depth output from top to bottom automatically satisfy the boundary conditions at the interlayer interfaces. The resulting horizontal electric field can be directly used for three-dimensional mesh boundary assignment without interpolation transformation. While achieving numerical stability throughout the process, the depth electric field profile required for three-dimensional magnetotelluric simulation is fully output.

[0139] Furthermore, based on the attenuation mode coefficients, the normalized propagation matrices of each interface, and the logarithmic scaling information, the horizontal electromagnetic field components at each interface are obtained as follows: For the target interface, the field state without scaling is calculated using the normalized propagation matrix, the base field matrix, and the attenuation mode coefficients corresponding to the interface. The relative scale factor is determined based on the difference between the cumulative logarithmic scale corresponding to the interface and the cumulative logarithmic scale corresponding to the surface interface. The field state without scaling is scaled using the relative scale factor to obtain the horizontal electromagnetic field components at the interface.

[0140] Specifically, when the substrate attenuation mode coefficients after the absorption scale are obtained... Then, using this coefficient, along with the normalized propagation matrix and logarithmic scale information of each interface, the horizontal electromagnetic field components at all interfaces are reconstructed layer by layer from top to bottom for subsequent generation of three-dimensional boundary conditions.

[0141] For any interface ,in, Read the normalized propagation matrix stored in the corresponding interface. , basis field matrix and the obtained attenuation mode coefficients Calculate the field state without scaling:

[0142] in, Let be a 4×1 complex vector containing the four field components at the interface that have not been scaled down. Since the value of is controlled to be around 1, the product result will not produce numerical overflow, but its value does not yet reflect the true physical field amplitude.

[0143] Then read the cumulative logarithmic scale corresponding to that interface from the storage. and the surface interface Corresponding cumulative logarithmic scale Calculate the difference between the two. And using this difference as an exponent, the relative scale factor is determined. .because It is the uppermost interface of the entire model, and physically, it has the longest propagation path from the base to the surface, resulting in the largest cumulative scale. Therefore, it usually has... This means the relative scale factor It is a positive number no greater than 1, which is completely within the safe range of floating-point representation and will not cause any overflow risk.

[0144] Finally, multiplying the aforementioned relative scale factor by the unscaled field state yields the true and complete horizontal electromagnetic field components at the interface:

[0145] in, The four components are as follows: , , and , representing the horizontal electric field in the x-direction, the horizontal electric field in the y-direction, the horizontal magnetic field in the x-direction, and the horizontal magnetic field in the y-direction at that depth, respectively. By changing , The values ​​are taken from top to bottom, starting from the surface. To the top of the base By performing the above calculations interface by interface, the horizontal electromagnetic field profiles at all N+1 layer interfaces can be obtained. This process only requires calculation of the relative scale. No need to restore the global scale factor that may overflow throughout the process. .

[0146] By using the normalized propagation matrix and the base field matrix corresponding to the target interface, along with the attenuation mode coefficients, the original field state without scale correction can be restored first. A relative scale factor is constructed by using the difference between the cumulative logarithms of the target interface and the surface interface. Relying solely on the difference in exponential calculations can avoid the over-limit problem caused by global large-scale numerical calculations. After combining the unscaled field state with the relative scale factor for correction, the complete horizontal electromagnetic field components of the interface can be directly obtained. The single-depth field solution is completed independently by relying on the logarithmic parameters stored independently for each interface. The output electromagnetic field naturally satisfies the interlayer continuity physical condition. The obtained horizontal electric field can be directly used for the layout of the three-dimensional mesh boundary without additional interpolation, conversion, or other post-processing steps, simplifying the calculation process.

[0147] In some embodiments, the solution method further includes: The horizontal electromagnetic field components at each interface are written into the top boundary edge, bottom boundary edge, and corresponding sidewall boundary edge of the three-dimensional staggered mesh according to their depth position to form a known boundary field vector. Based on the known boundary field vector, the preset three-dimensional internal field solution equation is solved to obtain the internal electric field under each polarization direction. The known boundary field vector is superimposed with the internal electric field to obtain the electric field distribution inside the three-dimensional computational domain.

[0148] After reconstructing the horizontal electromagnetic field components at all layer interfaces, these one-dimensional field profiles are applied to three-dimensional magnetotelluric (MT) forward modeling (MT) to provide physically consistent boundary conditions for the three-dimensional computational domain and drive the solution of the internal fields. Three-dimensional MT forward modeling typically applies physically valid boundary conditions to the outer boundary of the computational domain. The boundaries of the three-dimensional computational domain include the top boundary, bottom boundary, and four sidewall boundaries. Due to computational resource limitations, it is not possible to perform equal-scale partitioning of the entire space; therefore, appropriate boundary conditions need to be applied at the truncated boundaries to eliminate spurious reflections. For the medium near the side boundaries, it can usually be approximated as a one-dimensional layered column with conductivity varying only with depth. The one-dimensional field profile calculated based on this can be directly used as the known electric field values ​​of the three-dimensional boundary.

[0149] For each horizontal grid location on the boundary of the three-dimensional computational domain The local one-dimensional conductivity column at a given location is extracted based on the area or volume weight of adjacent cells. Specifically, for several three-dimensional conductivity tensor sampling points located around the boundary edge in the three-dimensional mesh, a weighted average is performed using the inverse of the distance to the boundary edge or the cell volume percentage as weights to obtain the one-dimensional conductivity distribution with depth at that boundary location. In one implementation, the conductivity components of the three-dimensional model are weighted and averaged according to the area around the boundary edge to construct an axis-aligned local one-dimensional column, and its off-diagonal components are set to zero. In another implementation, if the input is a full tensor containing vertically coupled components, vertical elimination is performed before entering the one-dimensional solver.

[0150] After obtaining the local one-dimensional column at each location on the boundary, the above steps are performed separately for the two linearly independent surface polarization directions, i.e., calculating the one-dimensional field profiles under the x-direction polarization and the y-direction polarization, respectively. The two linearly independent surface polarizations are: ,

[0151] The horizontal electric field components calculated for each local one-dimensional column under the two polarization directions. According to its depth position and horizontal grid position Write them separately into the staggered mesh edge assembly Specifically, it includes: Write at the top boundary edge The horizontal electric field component at that point is used as the boundary value of the incident field at the Earth's surface; at the bottom boundary edge, it is written... The horizontal electric field component at the location is used as the boundary value propagating upwards from the top of the substrate; at the edge of the sidewall boundary, the value is determined according to the depth layer at which the boundary edge is located. Write to the corresponding depth The horizontal electric field component at that location. Simultaneously, for edges adjacent to the boundary within the 3D mesh, the electric field value at the corresponding depth can be written in through interpolation or direct assignment, depending on the required computational accuracy. This forms the known boundary field vector for each polarization direction. This vector stores the known electric field value on the outer boundary edge and sets it to zero on the remaining internal degrees of freedom.

[0152] Then, the internal fields are solved using the same discrete operators as those used in the 3D internal forward modeling. Specifically, a time factor is employed. Conventions and discrete formats are used to obtain three-dimensional discrete electric field operators. This operator is obtained by discretizing the frequency domain electric field equations using the finite element method or the finite difference method. Given the boundary field vector... Applying this discrete operator, calculate The result serves as the equivalent excitation source driving the solution of the internal field, constituting the right-hand side of the three-dimensional internal field solution equation.

[0153] Solve the following three-dimensional internal field equations:

[0154] Obtain the internal unknown fields under each polarization direction The internal field is zero at the outer boundary and is uniquely determined by the equivalent excitation source in its internal degrees of freedom. Superimposing the known boundary field vector with the internal electric field:

[0155] Obtain the complete electric field distribution within the three-dimensional computational domain. The total field satisfies the given polarization boundary conditions on the boundary and the frequency domain electric field equation inside, and is the final output of the three-dimensional MT forward modeling.

[0156] By filling the corresponding horizontal electromagnetic fields at various depths onto the edges of various boundaries in a three-dimensional staggered grid and assembling the boundary field vectors, the boundary constraint format of the three-dimensional magnetotelluric discrete grid can be directly matched without additional conversion or interpolation processing. Based on this boundary field, the excitation terms of the three-dimensional equations can be constructed to solve for the unknown internal electric field, achieving physical self-consistency between the boundary constraints and the global field control equations. Finally, by superimposing the boundary field and the solved internal electric field, the complete global electric field distribution can be obtained, thus realizing the entire process from one-dimensional stable field calculation to three-dimensional forward simulation. The output electric field data can be directly used for the extraction of electrical parameters such as apparent resistivity and impedance, as well as for geological structure interpretation.

[0157] Furthermore, after completing the three-dimensional global electric field calculation, COMMEMI can also be used. 3D 2. The results were verified using a standard benchmark model. The apparent resistivity and impedance phase were compared to verify the overall effectiveness of the algorithm. In the low-frequency and mid-frequency range, the one-dimensional boundary field driven three-dimensional solution of this application can achieve high matching accuracy. The deviation in the high-frequency scenario comes from the three-dimensional numerical solution process such as mesh discretization and interpolation, and does not mean that the one-dimensional solution algorithm of this application has defects.

[0158] For example, refer to Figures 3-10 ,along The profile was constructed with 54 surface receivers, and four characteristic frequencies were selected for analysis, namely: and It covers a wide frequency range from low to high frequencies. The comparison metric is the two off-diagonal impedance components. Corresponding apparent resistivity (unit: Ω) (m) and impedance phase (unit: °). The comparison adopts the published results of EM3DANI three-dimensional anisotropic forward modeling by Han et al. in 2018.

[0159] For the apparent resistivity comparison at each frequency, the median absolute relative error at the receiving point is defined as follows: First, the absolute value of the relative error between the simulated value and the reference value at each receiving point is calculated, and then the median value of the relative errors at all 54 receiving points is taken. For impedance phase, since phase itself is an angular quantity, the absolute value of the absolute error between the simulated value and the reference value at each receiving point is defined, and then the median value is taken.

[0160] Reference Figures 3-8 In the low-frequency to mid-frequency range ( ), The median absolute relative error of apparent resistivity for both off-diagonal impedance components is less than approximately 1.6%, and the median absolute error of impedance phase is less than 0.5°. The one-dimensional boundary electric field generated by this invention can accurately drive three-dimensional simulations. The simulation results highly coincide with the standard reference solution across the entire profile, and there are no systematic deviations caused by numerical overflow in the one-dimensional boundary conditions. This accuracy fully meets the engineering requirements of practical geophysical exploration.

[0161] Reference Figure 9 and Figure 10At a high frequency of 100Hz, the median apparent resistivity difference is approximately 17%, and the median phase difference is approximately 4.8°. The main sources of high-frequency error include: insufficient discretization accuracy of the near-surface 3D mesh; the shallow skin depth of high-frequency electromagnetic waves, making them more sensitive to near-surface mesh size; the potential introduction of local errors through the field value interpolation method at the surface receiver points; and the use of area-weighted averaging for the one-dimensional cylinders on the boundary, which may not fully reflect the 3D effect in regions with drastic lateral changes in conductivity. These deviations all originate from the 3D numerical discretization and boundary simplification processes and are not inherent accuracy defects in the one-dimensional stable solution algorithm of this application. High-frequency error can be further reduced by refining the 3D mesh and improving the construction of local one-dimensional cylinders.

[0162] Thus, the solution method of this application possesses sufficient engineering accuracy within the core frequency band of low- to mid-frequency exploration, reliably supporting three-dimensional magnetotelluric numerical simulations. High-frequency errors can be further reduced by refining the three-dimensional mesh and detailing local one-dimensional columns. It avoids a second numerical overflow when multiple finite-layer matrices are multiplied consecutively, and can output the horizontal electric field at each depth after a given surface polarization, not just the surface impedance, thus directly serving the generation of three-dimensional magnetotelluric boundary conditions.

[0163] Although embodiments of this application have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting this application. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of this application.

Claims

1. A method for solving the stability of the magnetotelluric field based on exponential scaling and renormalization, characterized in that, It includes: Obtain the three-dimensional conductivity tensor and thickness of each layer in the layered medium model, and construct the two-dimensional equivalent horizontal conductivity tensor of each layer based on the three-dimensional conductivity tensor of each layer. Based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor, the exponential growth characteristic of the corresponding layer is determined, and the original field propagation matrix of the corresponding layer is scaled based on the exponential growth characteristic to generate the stabilization layer matrix of each layer. Starting from the bottom layer of the layered medium model, the intermediate matrix corresponding to the current layer is obtained by going up layer by layer based on the stabilization layer matrix of the current layer and the propagation matrix accumulated by the lower layer. The intermediate matrix corresponding to the current layer is normalized and the corresponding logarithmic scale information is recorded to obtain the normalized propagation matrix of each layer interface and its corresponding cumulative logarithmic scale. Based on the given surface electric field boundary conditions and the field attenuation conditions of the base half-space, the attenuation mode coefficients in the base half-space are obtained using the normalized propagation matrix of the topmost interface and the logarithmic scale information. Based on the attenuation mode coefficients, the normalized propagation matrices of each interface, and the logarithmic scaling information, the horizontal electromagnetic field components at each interface are obtained.

2. The solution method according to claim 1, characterized in that, The determination of the exponential growth characteristic quantities of different electromagnetic modes in the corresponding layer based on the thickness of each layer and the two-dimensional equivalent horizontal conductivity tensor includes: The two-dimensional equivalent horizontal conductivity tensor of the corresponding layer is subjected to eigenvalue decomposition to obtain the principal conductivity eigenvalues ​​corresponding to the two modes of the corresponding layer, wherein the two modes include the first mode and the second mode. Based on the two principal conductivity eigenvalues, determine the complex vertical propagation constants corresponding to the two principal conductivity eigenvalues ​​of the corresponding layer; Multiply the two complex vertical propagation constants by the thickness of the corresponding layer and take the real part to obtain the corresponding exponential growth characteristic.

3. The solution method according to claim 2, characterized in that, The scaling process of the original field propagation matrix of the corresponding layer based on the exponentially growing feature quantity to generate the stabilization layer matrix of each layer includes: Determine the original hyperbolic cosine function value and the original hyperbolic sine function value corresponding to each mode respectively; Using the exponentially growing feature, the original hyperbolic cosine function value and the original hyperbolic sine function value of each mode are exponentially scaled to obtain the scaled hyperbolic cosine function value and the scaled hyperbolic sine function value. Based on the complex vertical propagation constant and the corresponding principal conductivity eigenvalue of each mode, the characteristic impedance of each mode is determined respectively; Based on the scaled hyperbolic cosine function value, the scaled hyperbolic sine function value, and the characteristic impedance of each mode, a scaling layer matrix of the corresponding layer in the local principal axis coordinate system is constructed. Transform the scaling layer matrix to the global coordinate system to obtain the stabilization layer matrix of the corresponding layer.

4. The solution method according to claim 3, characterized in that, The scaling layer matrix of the corresponding layer in the local principal axis coordinate system is constructed based on the scaled hyperbolic cosine function value, the scaled hyperbolic sine function value, and the characteristic impedance of each mode, including: Based on the scaled hyperbolic cosine function value corresponding to the first mode, the scaled hyperbolic sine function value corresponding to the first mode, the scaled hyperbolic cosine function value corresponding to the second mode, and the scaled hyperbolic sine function value corresponding to the second mode, as well as the characteristic impedance and the inverse of the characteristic impedance, two sets of two-dimensional field transfer sub-matrices are constructed. The two sets of two-dimensional field transfer sub-matrices are merged to obtain the scaling layer matrix of the corresponding layer in the local principal axis coordinate system.

5. The solution method according to claim 3, characterized in that, The step of exponentially scaling the original hyperbolic cosine function values ​​and the original hyperbolic sine function values ​​of each mode using the exponentially growing feature to obtain scaled hyperbolic cosine function values ​​and scaled hyperbolic sine function values ​​includes: Compare the magnitudes of the two exponentially growing feature quantities in the corresponding layer, and select the maximum value as the common offset of the corresponding layer; For each mode in each layer, the complex vertical propagation constant of that mode is multiplied by the thickness of that layer to obtain the complex propagation thickness value of that mode; Construct positive exponent terms with the complex propagation thickness value as the exponent, and negative exponent terms with the negative value of the complex propagation thickness value as the exponent, respectively; Subtract the common offset from the exponent portion of the positive exponent term and the exponent portion of the negative exponent term respectively to obtain the scaled positive exponent term and the scaled negative exponent term. The hyperbolic cosine function value of this mode is obtained based on the sum of the scaled positive exponent term and the scaled negative exponent term. The modally scaled hyperbolic sine function value is obtained based on the difference between the scaled positive exponent term and the scaled negative exponent term.

6. The solution method according to claim 1, characterized in that, The intermediate matrix corresponding to the current layer is normalized, and the corresponding logarithmic scaling information is recorded, including: Extract the maximum value among all elements in the intermediate matrix and use it as the normalization factor for this normalization process; Divide each element of the intermediate matrix by the normalization factor to obtain the normalized propagation matrix of the current layer interface; The natural logarithm of the normalization factor and the common offset corresponding to the current layer are added to the cumulative logarithmic scale accumulated in the lower layer to obtain the logarithmic scale information of the current layer interface.

7. The solution method according to claim 1, characterized in that, The attenuation mode coefficients in the base half-space are obtained by using the normalized propagation matrix of the topmost interface and the logarithmic scaling information, based on the given surface electric field boundary conditions and the field attenuation conditions of the base half-space: Based on the electrical parameters of the base half-space, construct the base field matrix that retains only the two downward decaying electromagnetic modes in the base half-space; Using the normalized propagation matrix of the topmost interface, the logarithmic scaling information of the topmost interface, and the base field matrix, a mapping matrix from the attenuation mode coefficients to the horizontal electric field of the ground surface is constructed. The exponential factor of the cumulative logarithmic scale of the topmost interface is incorporated into the decay mode coefficients to be solved. The attenuation mode coefficients are obtained by solving the linear equations formed by the mapping matrix based on the given surface electric field boundary conditions.

8. The solution method according to claim 1, characterized in that, The horizontal electromagnetic field components at each interface, obtained based on the attenuation mode coefficients, the normalized propagation matrices of each interface, and the logarithmic scaling information, include: For the target interface, the field state without scaling is calculated using the normalized propagation matrix corresponding to the interface, the basis field matrix, and the attenuation mode coefficients. The relative scale factor is determined based on the difference between the cumulative logarithmic scale corresponding to the interface and the cumulative logarithmic scale corresponding to the surface interface. The field state before scaling is scaled using the relative scale factor to obtain the horizontal electromagnetic field component at the interface.

9. The solution method according to any one of claims 1-8, characterized in that, The construction of the two-dimensional equivalent horizontal conductivity tensor for each layer based on the three-dimensional conductivity tensor of each layer includes: The vertical conductivity component, horizontal conductivity component, and coupled conductivity component between the vertical and horizontal components are determined based on the three-dimensional conductivity tensor. Based on the characteristic that the electromagnetic field in a one-dimensional layered medium remains unchanged along the horizontal direction, and using the physical constraint that the vertical current is zero, the correlation between the vertical electric field and the horizontal electric field is established. Based on the aforementioned correlation, the vertical electric field component of the three-dimensional conductivity tensor is eliminated to obtain a two-dimensional equivalent horizontal conductivity tensor containing only the horizontal component. Each component of the two-dimensional equivalent horizontal conductivity tensor is obtained by subtracting a correction term from the horizontal conductivity component. The correction term is the quotient obtained by dividing the product of the coupled conductivity components by the vertical conductivity component.

10. The solution method according to any one of claims 1-8, characterized in that, The solution method also includes: The horizontal electromagnetic field components at each interface are written into the top boundary edge, bottom boundary edge, and corresponding sidewall boundary edge of the three-dimensional staggered grid according to the depth position to form a known boundary field vector. Based on the known boundary field vector, the preset three-dimensional internal field equation is solved to obtain the internal electric field under each polarization direction; The known boundary field vector is superimposed with the internal electric field to obtain the electric field distribution inside the three-dimensional computational domain.

Images