A method and system for measuring apparent resistivity of horizontally layered earth

Abstract

The application relates to the technical field and discloses a method and system for measuring and calculating apparent resistivity of horizontally layered earth, wherein a characteristic equation of horizontally layered earth is constructed according to an integral variable set obtained through an exponential discrete sampling method and horizontally layered earth parameters; an uplink electromagnetic wave amplitude kernel function of the characteristic equation at each integral variable is calculated to obtain an integral upper limit of horizontally layered earth, and the integral upper limit of horizontally layered earth is set as the maximum value of the integral variable satisfying a preset criterion; based on the integral upper limit of horizontally layered earth, the integral range of the constructed apparent resistivity relationship formula of an unbounded domain is limited to obtain an apparent resistivity relationship formula of a bounded domain; and based on a 151-point Gauss-Legendre integral algorithm and Euler transformation, the apparent resistivity relationship formula of the bounded domain is iteratively solved to obtain a measurement and calculation result of the apparent resistivity. The method meets the dual requirements of precision and efficiency for the measurement and calculation of the apparent resistivity of horizontally layered earth.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geological exploration, and in particular to a method and system for calculating apparent resistivity of horizontally layered earth. BACKGROUND

[0002] In the field of geological exploration, accurate calculation of the apparent resistivity of horizontally layered earth is of great significance for analyzing underground geological structure and exploring resources. The apparent resistivity of horizontally layered earth is the resistivity of each layer calculated according to the uniform earth resistivity measurement method when the electrical properties of underground media are uneven and the surface is uneven. It is influenced by many factors such as the true resistivity of underground media, electrode spacing, electrode arrangement, and measurement depth.

[0003] The existing methods for calculating the apparent resistivity of horizontally layered earth include the complex mirror image method and the 141 linear filtering method. The complex mirror image method needs to use the sampling law to obtain the dynamic waveform range of the first layer uplink electromagnetic wave amplitude kernel function with respect to the integral variable. There is no related method description at present. The complex mirror image method needs to use complex matrix calculation methods such as singular decomposition and projection transformation to determine the strength and mode of the complex mirror image, which has large calculation amount and low efficiency. The 141 linear filtering method uses logarithmic equidistant sampling method. In order to obtain the optimal filtering coefficient and exponential sampling coefficient, a large amount of calculation is needed, and when the working condition changes, the specific working condition needs to be recalculated, so the efficiency of calculation is low. Therefore, the existing method for calculating the apparent resistivity of horizontally layered earth cannot meet the dual requirements of precision and efficiency of the apparent resistivity calculation of horizontally layered earth.

[0004] Therefore, how to meet the dual requirements of precision and efficiency of the apparent resistivity calculation of horizontally layered earth has become a technical problem to be solved by those skilled in the art. SUMMARY

[0005] The present application provides a method and system for calculating the apparent resistivity of horizontally layered earth to solve the technical problem of how to meet the dual requirements of precision and efficiency of the apparent resistivity calculation of horizontally layered earth, and achieve the effect of meeting the dual requirements of precision and efficiency of the apparent resistivity calculation of horizontally layered earth.

[0006] In a first aspect, the present application provides a method for calculating the apparent resistivity of horizontally layered earth, which comprises:

[0007] The integral variable set is obtained by the exponential discrete sampling method, and the characteristic equation of horizontally layered earth is constructed according to the integral variable set and the obtained horizontal layered earth parameters of the target earth region.

[0008] According to the first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation under each integral variable, a horizontal layered earth integral upper limit is obtained, the horizontal layered earth integral upper limit is set as the maximum value of the integral variable meeting a preset criterion, the preset criterion is used to reflect the relationship between the local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set and the allowable error, and the allowable error is set as the fluctuation range determined according to the first layered uplink electromagnetic wave amplitude kernel function changing with the integral variable;

[0009] According to the obtained pole distance of the Wenner quadrupole method, an unbounded domain apparent resistivity relationship of the horizontal layered earth is constructed;

[0010] Based on the horizontal layered earth integral upper limit, the integral range of the unbounded domain apparent resistivity relationship is limited to obtain a bounded domain apparent resistivity relationship;

[0011] Based on the 151-point Gauss-Legendre integral algorithm and Euler transformation, the bounded domain apparent resistivity relationship is iteratively solved to obtain the apparent resistivity calculation result.

[0012] Preferably, the characteristic equation of the horizontal layered earth is constructed according to the integral variable set and the obtained horizontal layered earth parameters of the target large area, comprising:

[0013] The horizontal layered earth parameters of the target large area are obtained, the horizontal layered earth parameters include the soil resistivity of each layer and the z-axis coordinate of the layer boundary of each layer, and the reflection coefficient of each layer boundary is constructed according to the soil resistivity;

[0014] Based on the electromagnetic characteristics of the horizontal layered earth, the interlayer relationship of the layer boundary of the horizontal layered earth is determined according to the reflection coefficient, the integral variable set and the z-axis coordinate of the layer boundary;

[0015] Based on the electromagnetic field boundary condition, the characteristic equation of the horizontal layered earth is constructed, the characteristic equation is composed of a characteristic equation matrix, a characteristic equation vector and a to-be-solved characteristic vector, the elements of the characteristic equation matrix represent the interlayer relationship, the elements of the characteristic equation vector represent the known electromagnetic field boundary excitation condition, and the to-be-solved characteristic vector includes the first layered uplink electromagnetic wave amplitude kernel function.

[0016] Preferably, the first layered uplink electromagnetic wave amplitude kernel function of the characteristic equation under each integral variable is calculated to obtain the horizontal layered earth integral upper limit, comprising:

[0017] According to the integral variable set, the characteristic equation is iteratively solved to obtain the to-be-solved characteristic vector corresponding to each integral variable in the integral variable set, and the first layered uplink electromagnetic wave amplitude kernel function in each to-be-solved characteristic vector is obtained;

[0018] According to each first layered uplink electromagnetic wave amplitude kernel function, a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function is obtained, and an allowable error is constructed according to the fluctuation range.

[0019] According to each first layered uplink electromagnetic wave amplitude kernel function, a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set is obtained.

[0020] According to the local absolute error and the allowable error, a preset criterion is constructed, and the preset criterion is solved to obtain a maximum value of the integral variable that satisfies the preset criterion, and the maximum value is taken as a horizontal layered earth integral upper limit.

[0021] Preferably, the 151-point Gauss-Legendre integral algorithm and Euler transformation are used to iteratively solve the bounded domain apparent resistivity relationship to obtain an apparent resistivity calculation result, including:

[0022] The 151-point Gauss-Legendre integral algorithm is used to solve the bounded domain apparent resistivity relationship to obtain an apparent resistivity Gauss calculation formula.

[0023] An Euler transformation coefficient function is introduced into the apparent resistivity Gauss calculation formula to obtain an apparent resistivity Euler transformation calculation formula.

[0024] The apparent resistivity Euler transformation calculation formula is iteratively solved to obtain an apparent resistivity calculation result.

[0025] Preferably, the apparent resistivity Euler transformation calculation formula is expressed as:

[0026]

[0027] wherein, denotes a pole distance, denotes a pole distance number of the Wenner quadrupole method, denotes a pole distance corresponding apparent resistivity, denotes a soil resistivity of a first layer, denotes an integral point index of the 151-point Gauss-Legendre integral, denotes a position coefficient of the 151-point Gauss-Legendre integral, denotes a horizontal layered earth integral upper limit, denotes a weight coefficient of 151-point Gauss-Legendre integration, denotes a coefficient function of Euler transformation, denotes an integral variable corresponding first layered uplink electromagnetic wave amplitude kernel function, denotes corresponding first type zero order Bessel function, denotes corresponding first type zero order Bessel function.

[0028] In a second aspect, the present application also provides a system for calculating apparent resistivity of horizontally layered earth, which implements the method for calculating apparent resistivity of horizontally layered earth as described above, and comprises a characteristic equation construction unit, a horizontally layered earth integral upper limit calculation unit, an apparent resistivity relationship construction unit, an apparent resistivity relationship conversion unit and an apparent resistivity calculation unit.

[0029] The characteristic equation construction unit is configured to obtain an integral variable set by an exponential discrete sampling method, and construct a characteristic equation of horizontally layered earth according to the integral variable set and acquired horizontally layered earth parameters of a target earth region.

[0030] The horizontally layered earth integral upper limit calculation unit is configured to calculate a first layered uplink electromagnetic wave amplitude kernel function of the characteristic equation under each integral variable, and obtain a horizontally layered earth integral upper limit, which is set as a maximum value of the integral variable satisfying a preset criterion, and the preset criterion is configured to reflect a relationship between a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to a highest integral variable in the integral variable set and an allowable error, which is set according to a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with respect to the integral variable.

[0031] The apparent resistivity relationship construction unit is configured to construct an unbounded domain apparent resistivity relationship of the horizontally layered earth according to an acquired pole distance of Wenner quadrupole method.

[0032] The apparent resistivity relationship conversion unit is configured to limit an integral range of the unbounded domain apparent resistivity relationship based on the horizontally layered earth integral upper limit, and obtain a bounded domain apparent resistivity relationship.

[0033] The apparent resistivity calculation unit is configured to iteratively solve the bounded domain apparent resistivity relationship based on a 151-point Gauss-Legendre integral algorithm and Euler transformation, and obtain an apparent resistivity calculation result.

[0034] Preferably, the characteristic equation construction unit comprises a reflection coefficient calculation module, an interlayer relationship construction module and a characteristic equation construction module.

[0035] The reflection coefficient calculation module is configured to acquire horizontal layered earth parameters of a target large area, the horizontal layered earth parameters comprising soil resistivity of each layer and z-axis coordinates of layer boundaries of each layer, and construct reflection coefficients of the layer boundaries according to the soil resistivity.

[0036] The interlayer relationship construction module is configured to determine an interlayer relationship of the layer boundaries of the horizontal layered earth based on electromagnetic characteristics of the horizontal layered earth, the reflection coefficients, the integral variable set and the z-axis coordinates of the layer boundaries.

[0037] The characteristic equation construction module is configured to construct a characteristic equation of the horizontal layered earth based on electromagnetic field boundary conditions, the characteristic equation comprising a characteristic equation matrix, a characteristic equation vector and a to-be-solved eigenvector, elements of the characteristic equation matrix representing the interlayer relationship, elements of the characteristic equation vector representing known electromagnetic field boundary excitation conditions, and the to-be-solved eigenvector comprising the first-layer uplink electromagnetic wave amplitude kernel function.

[0038] Preferably, the horizontal layered earth integral upper limit calculation unit comprises a first calculation module, an allowable error construction module, a local absolute error construction module and a second calculation module.

[0039] The first calculation module is configured to iteratively solve the characteristic equation according to the integral variable set to obtain the to-be-solved eigenvector corresponding to each integral variable in the integral variable set, and acquire the first-layer uplink electromagnetic wave amplitude kernel function in each to-be-solved eigenvector.

[0040] The allowable error construction module is configured to obtain a fluctuation range of the first-layer uplink electromagnetic wave amplitude kernel function according to each first-layer uplink electromagnetic wave amplitude kernel function, and construct an allowable error according to the fluctuation range.

[0041] The local absolute error construction module is configured to obtain a local absolute error of the first-layer uplink electromagnetic wave amplitude kernel function corresponding to a highest integral variable in the integral variable set according to each first-layer uplink electromagnetic wave amplitude kernel function.

[0042] The second calculation module is configured to construct a preset criterion according to the local absolute error and the allowable error, and solve the preset criterion to obtain a maximum value of the integral variable satisfying the preset criterion, and take the maximum value as a horizontal layered earth integral upper limit.

[0043] Preferably, the apparent resistivity calculation unit comprises a third calculation module, an Euler transformation module and a fourth calculation module.

[0044] The third calculation module is configured to solve the bounded domain apparent resistivity relationship by using a 151-point Gauss-Legendre integral algorithm to obtain an apparent resistivity Gauss calculation formula.

[0045] The Euler transformation module is configured to introduce an Euler transformation coefficient function into the apparent resistivity Gauss calculation formula to obtain an apparent resistivity Euler transformation calculation formula.

[0046] The fourth calculation module is configured to iteratively solve the apparent resistivity Euler transformation calculation formula to obtain an apparent resistivity calculation result.

[0047] Preferably, the apparent resistivity Euler transformation calculation formula is expressed as:

[0048]

[0049] wherein, denotes a pole distance, denotes a pole distance number of the Wenner quadrupole method, denotes a pole distance corresponding apparent resistivity, denotes a soil resistivity of a first layer, denotes an integral point index of the 151-point Gauss-Legendre integral, denotes a position coefficient of the 151-point Gauss-Legendre integral, denotes an upper limit of horizontal layered earth integral, denotes a weight coefficient of the 151-point Gauss-Legendre integral, denotes a coefficient function of Euler transformation, denotes an integral variable corresponding uplink electromagnetic wave amplitude kernel function of the first layer, denotes corresponding first type zero-order Bessel function, denotes corresponding first type zero-order Bessel function.

[0050] The present application provides a horizontal layered earth apparent resistivity calculation method and system, and the beneficial effects of the present application are as follows:

[0051] The apparent resistivity measurement method of the horizontally layered earth disclosed in the application adopts the equal exponent sampling method to determine the integral variable, can obtain the sharp change of the first layered uplink electromagnetic wave amplitude kernel function in the low frequency band, and can obtain the asymptotic value of the first layered uplink electromagnetic wave amplitude kernel function in the high frequency band, so as to obtain the dynamic change waveform of the first layered uplink electromagnetic wave amplitude kernel function, and ensure the accuracy of the measurement. Meanwhile, the exponent discretization step of the integral variable is determined, avoiding the influence of the soil parameter change. The upper limit of the horizontally layered earth integral is determined by the soil parameter, reflecting the difference of the soil parameter, and the optimal solution solving of the filtering coefficient and the exponent sampling coefficient is not needed by using a large number of cases, which significantly reduces the complexity of the measurement. The Euler transformation of the apparent resistivity is realized by using the 151-point Gauss-Legendre numerical integral algorithm combined with the Euler transformation coefficient, which can effectively improve the measurement efficiency of the apparent resistivity of the horizontally layered earth and reduce the complexity of the measurement. BRIEF DESCRIPTION OF DRAWINGS

[0052] Figure 1 is a horizontal layered earth apparent resistivity measurement method step schematic diagram provided by one preferred embodiment of the application;

[0053] Figure 2 is a horizontally layered earth model diagram provided by one preferred embodiment of the application;

[0054] Figure 3 is a structure schematic diagram of a horizontally layered earth apparent resistivity measurement system provided by one preferred embodiment of the application;

[0055] REFERENCE NUMERALS:

[0056] 1-feature equation construction unit, 2-horizontally layered earth integral upper limit calculation unit, 3-apparent resistivity relationship construction unit, 4-apparent resistivity relationship conversion unit, 5-apparent resistivity calculation unit. DETAILED DESCRIPTION

[0057] The embodiments of the application will be described in detail below with reference to the drawings. The embodiments are given only for the purpose of illustration and should not be understood as limiting the application. The accompanying drawings are used for reference and illustration only and do not constitute a limitation on the scope of patent protection. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application. In the description of the application, the terms "first", "second", "third" and the like are used only for the purpose of description and should not be understood as indicating or implying relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features limited by "first", "second", "third" and the like can explicitly or implicitly include one or more features. In the description of the application, unless otherwise specified, the meaning of "multiple" is two or more.

[0058] In the description of the present application, it should be noted that, unless otherwise explicitly defined and limited, the terms "mounting", "connecting", "connecting" should be understood broadly, for example, it can be fixedly connected, or it can be detachably connected, or integrally connected, it can be mechanically connected, or it can be electrically connected, it can be directly connected, or indirectly connected through an intermediate medium, it can be the communication inside two elements. The terms "vertical", "horizontal", "left", "right", "up", "down" and similar expressions used herein are for illustrative purposes only, and do not indicate or imply that the device or element referred to must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be construed as limiting the present application. The term "and / or" used herein includes any and all combinations of one or more related listed items. The specific meaning of the above terms in the present application can be understood in specific cases by those of ordinary skill in the art.

[0059] In the description of the present application, it should be noted that, unless otherwise defined, all technical and scientific terms used in the present application have the same meaning as understood by those skilled in the art. The terms used in the specification of the present application are only for the purpose of describing the specific embodiments, and are not intended to limit the present application. The specific meaning of the above terms in the present application can be understood in specific cases by those of ordinary skill in the art.

[0060] Please refer to Figure 1 The method for measuring and calculating the apparent resistivity of horizontally layered earth is shown in the step schematic diagram of the method for measuring and calculating the apparent resistivity of horizontally layered earth, and in the embodiments of the present application, a method for measuring and calculating the apparent resistivity of horizontally layered earth is provided, which comprises:

[0061] S1, obtaining an integral variable set by an exponential discrete sampling method, and constructing a characteristic equation of horizontally layered earth according to the integral variable set and the horizontal layered earth parameters of the target large area region obtained; the technical scheme of the present application is implemented under the condition that the horizontal layered earth parameters and the four-pole method pole distance arrangement are known, and the horizontal layered earth parameters include the soil resistivity of each layer, the thickness of each layer and the boundary of each layer.

[0062] Further, the integral variable is determined by using the equal exponential sampling method, the integral variable is set as the integral variable, the integral variable forms an integral variable set, and the integral variable set is represented as follows:

[0063]

[0064] Wherein, represents the integral variable, represents the integral variable exponential discretization step index.

[0065] Based on the electromagnetic characteristics of the horizontally layered ground, the boundary conditions of electromagnetic wave propagation in the layered medium of the horizontally layered ground are converted into matrix equations by physical modeling and mathematical derivation to construct a characteristic equation. To describe the electromagnetic wave propagation characteristics of each layer of the horizontally layered ground, a Hankel integral kernel function is used to represent the amplitude characteristics of the uplink electromagnetic wave in the layer, which is an uplink electromagnetic wave amplitude kernel function, and a Hankel integral kernel function is used to represent the amplitude characteristics of the downlink electromagnetic wave in the layer, which is a downlink electromagnetic wave amplitude kernel function.

[0066] Further, according to the soil resistivity, the reflection coefficient of each layered boundary is constructed, and the reflection coefficient of each layered boundary is represented as:

[0067]

[0068] wherein, represents a layered index, represents the reflection coefficient of the first layered, represents the soil resistivity of the first layered, represents the soil resistivity of the first layered, represents the soil resistivity of the first layered. The upper and lower boundaries of each layer of the horizontally layered ground satisfy the electromagnetic field boundary conditions, the uplink electromagnetic wave amplitude kernel function and the uplink electromagnetic wave amplitude kernel function of the adjacent layer satisfy a linear relationship, the linear relationship is represented, and a characteristic equation is obtained, and the characteristic equation is represented as:

[0069]

[0070]

[0071] wherein, represents a characteristic equation matrix composed of the interlayer relationship of the layered boundary of the horizontally layered ground, represents a to-be-solved characteristic vector composed of the uplink electromagnetic wave amplitude kernel function and the downlink electromagnetic wave amplitude kernel function of the horizontally layered ground, represents a characteristic equation vector composed of known electromagnetic field boundary excitation conditions.

[0072] The to-be-solved characteristic vector is represented as:

[0073]

[0074] wherein, represents the first layered downlink electromagnetic wave amplitude kernel function, represents the first layered uplink electromagnetic wave amplitude kernel function, represents the number of layers, represents the downlink electromagnetic wave amplitude kernel function of the first layer, represents the downlink electromagnetic wave amplitude kernel function of the first layer, represents the downlink electromagnetic wave amplitude kernel function of the first layer.​​ Layer upgoing electromagnetic wave amplitude kernel function.

[0075] Characteristic equation matrix is is a matrix of dimension n x n, whose elements are given by:

[0076]

[0077]

[0078]

[0079]

[0080]

[0081]

[0082] where, is the coefficient of the surface normal boundary condition with respect to is the z-axis coordinate of the boundary of the i-th layer, is the coefficient of the surface normal boundary condition with respect to is the coefficient of the convergence boundary condition of the last layer of the layered earth with respect to is the coefficient of the surface normal boundary condition with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the potential boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to is the coefficient of the normal current density continuity boundary condition of the i-th layer of the layered earth with respect to

[0083] Characteristic equation vector is an n x 1 column vector whose non-zero elements are:

[0084]

[0085] where, ​​denotes the column vector to the first element, denotes the column vector to the second element.

[0086] S2, according to the first layered uplink electromagnetic wave amplitude kernel function corresponding to the feature equation under each integral variable, a horizontal layered earth integral upper limit is obtained, the horizontal layered earth integral upper limit is set as the maximum value of the integral variable which satisfies the preset criterion, the preset criterion is used to reflect the relationship between the local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set and the allowable error, the allowable error is set as determined according to the fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with the integral variable; in the preferred embodiment of the present application, the first layered uplink electromagnetic wave amplitude kernel function of the feature equation of the horizontal layered earth integral variable is calculated by iteration, and the preset criterion judgment is performed on the first layered uplink electromagnetic wave amplitude kernel function, so as to determine the horizontal layered earth integral upper limit.

[0087] Let , the feature equation is solved, and the integral variable is The corresponding to-be-solved characteristic vector is obtained, and is recorded.

[0088] Let , the feature equation is solved iteratively until The to-be-solved characteristic vector corresponding to each integral variable in the integral variable set is obtained, and corresponding to the to-be-solved characteristic vector is recorded.

[0089] According to each first layered uplink electromagnetic wave amplitude kernel function, the fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with the integral variable is obtained, and is expressed as follows:

[0090]

[0091] , wherein, denotes the fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function, denotes the maximum value of the first layered uplink electromagnetic wave amplitude kernel function, denotes the minimum value of the first layered uplink electromagnetic wave amplitude kernel function.

[0092] Further, according to the fluctuation range, the allowable error is constructed, according to each first layered uplink electromagnetic wave amplitude kernel function, the local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set is obtained, and according to the allowable error and the local absolute error, the preset criterion is constructed, and the preset criterion is expressed as follows:

[0093]

[0094] wherein, represents the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable, represents the integral variable satisfying the preset criterion, represents the integral variable index corresponding to the integral variable satisfying the preset criterion.

[0095] According to the preset criterion, the maximum value of the integral variable satisfying the preset criterion is obtained, and the maximum value is taken as the horizontal layered earth integral upper limit.

[0096] In the preferred embodiment of the present application, the equal index sampling method is used to determine the integral variable, which can obtain the sharp change of the first layered uplink electromagnetic wave amplitude kernel function in the low frequency band, and can also obtain the asymptotic value of the first layered uplink electromagnetic wave amplitude kernel function in the high frequency band, so as to obtain the dynamic change waveform of the entire first layered uplink electromagnetic wave amplitude kernel function, and ensure the accuracy of the measurement. At the same time, the integral variable index discretization step is determined, avoiding the influence of soil parameter change. The horizontal layered earth integral upper limit is determined by the soil parameter, which reflects the difference of the soil parameter, and does not need to use a large number of cases to solve the optimal solution of the filtering coefficient and the index sampling coefficient, which significantly reduces the complexity of the measurement.

[0097] S3, according to the obtained pole distance of the Wenner four-pole method, the apparent resistivity relationship formula of the unbounded domain of the horizontal layered earth is constructed; in the preferred embodiment of the present application, the measurement of the horizontal layered earth parameter adopts the Wenner four-pole method, including 2 power supply electrodes and 2 measurement electrodes, which are symmetrically arranged between the four electrodes. The distance between the electrodes is the pole distance, and a total of measurements are performed, corresponding to groups of pole distances.

[0098] The apparent resistivity calculation formula of the pole distance is constructed, and the apparent resistivity calculation formula is represented as:

[0099]

[0100] wherein, represents the pole distance, represents the pole distance corresponding to the apparent resistivity, represents the potential function.

[0101] As Figure 2 shown in the horizontal layered earth model diagram, the potential function is analyzed according to the horizontal layered earth model diagram to obtain the potential function analytical expression:

[0102]

[0103]

[0104] wherein, denotes the corresponding first kind zero order Bessel function, denotes the potential function integral variable, denotes the first layered soil resistivity, denotes the electrode spacing the corresponding first kind zero order Bessel function.

[0105] Since the range of the potential function integral variable of the potential function analytical expression is an infinite bound range, the apparent resistivity obtained according to the potential function analytical expression is an unbounded domain apparent resistivity relationship, and the unbounded domain apparent resistivity relationship is:

[0106]

[0107] S3, limiting the integral range of the unbounded domain apparent resistivity relationship based on the horizontal layered geodetic integral upper limit, to obtain a bounded domain apparent resistivity relationship; the upper limit of the potential function integral variable of the unbounded domain apparent resistivity relationship is set to the horizontal layered geodetic integral upper limit, to obtain a bounded domain apparent resistivity relationship, and the bounded domain apparent resistivity relationship is expressed as:

[0108]

[0109] wherein, denotes the horizontal layered geodetic integral upper limit.

[0110] S4, based on the 151-point Gauss-Legendre integral algorithm and Euler transformation, the bounded domain apparent resistivity relationship is iteratively solved to obtain the apparent resistivity calculation result; in the preferred embodiment of the present application, the 151-point Gauss-Legendre integral algorithm and Euler transformation are combined to iteratively solve the bounded domain apparent resistivity relationship, the 151-point Gauss-Legendre integral algorithm is a specific implementation form of the Gauss-Legendre integral algorithm using 151 nodes, and the 151-point Gauss-Legendre integral algorithm is used to solve the bounded domain apparent resistivity relationship to obtain an apparent resistivity Gauss calculation formula, and the apparent resistivity Gauss calculation formula is:

[0111]

[0112] wherein, denotes the integral point index of the 151-point Gauss-Legendre integral algorithm, denotes the corresponding first kind zero order Bessel function, denotes the corresponding first kind zero order Bessel function, denotes the weight coefficient of the 151-point Gauss-Legendre integral algorithm, The specific values of the position coefficient and the weight coefficient of the 151-point Gauss-Legendre integral algorithm need to be determined by the position coefficient and weight coefficient table of the 151-point Gauss-Legendre integral algorithm, which is public content and will not be specifically shown in the present application.

[0113] The Euler transformation coefficient function is introduced into the apparent resistivity calculation formula to obtain an apparent resistivity Euler transformation calculation formula, and the apparent resistivity Euler transformation calculation formula is:

[0114]

[0115] wherein, The coefficient function of the Euler transformation is represented as f(x), and then

[0116]

[0117] wherein, The complementary error function is represented as erfc(x), The integral variable is represented as x, The corresponding first layered uplink electromagnetic wave amplitude kernel function is represented as The first layered uplink electromagnetic wave amplitude kernel function is represented as The first type zero-order Bessel function is represented as J0(x), The first type zero-order Bessel function is represented as J0(x). The first type zero-order Bessel function is represented as J0(x).

[0118] For the in the apparent resistivity Euler transformation calculation formula The method for solving the first layered uplink electromagnetic wave amplitude kernel function in S1 is still adopted, and the difference lies in that the integral variable The value of x is , so the specific solving process will not be described again.

[0119] Let the pole distance number be , and the apparent resistivity Euler transformation calculation formula is solved to obtain the apparent resistivity measurement result of the first four-pole method.

[0120] Let , and the apparent resistivity Euler transformation calculation formula is iteratively solved until > , to obtain the final apparent resistivity measurement result.

[0121] In the preferred embodiments of the present application, the 151-point Gauss-Legendre numerical integral algorithm is used in combination with the Euler transformation coefficient to realize the Euler transformation of the apparent resistivity, which can effectively improve the measurement efficiency of the apparent resistivity of the horizontal layered earth and reduce the complexity of the measurement.

[0122] The apparent resistivity calculation method of the horizontal layered earth of the application is compared with the calculation results of the complex mirror image method and the 141-point linear filtering method, and a 10-layer horizontal layered earth example is verified. The horizontal layered earth parameters of the 10-layer horizontal layered earth are as follows: the soil resistivity of the first layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the first layer and the soil of the second layer is 1 m; the soil resistivity of the second layer is 1000 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the second layer and the soil of the third layer is 2 m; the soil resistivity of the third layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the third layer and the soil of the fourth layer is 3 m; the soil resistivity of the fourth layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the fourth layer and the soil of the fifth layer is 4 m; the soil resistivity of the fifth layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the fifth layer and the soil of the sixth layer is 5 m; the soil resistivity of the sixth layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the sixth layer and the soil of the seventh layer is 6 m; the soil resistivity of the seventh layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the seventh layer and the soil of the eighth layer is 5 m; the soil resistivity of the eighth layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the eighth layer and the soil of the ninth layer is 5 m; the soil resistivity of the ninth layer is 100 Ω·m, the thickness is 1 m, the boundary z-axis coordinate of the soil of the ninth layer and the soil of the tenth layer is 9 m; the soil resistivity of the tenth layer is 100 Ω·m, the thickness is 1 m.

[0123] The calculation results of the complex mirror image method, the 141-point linear filtering method and the apparent resistivity calculation method of the horizontal layered earth of the application are shown in the following table:

[0124]

[0125] As can be seen from Table 1, the error of the calculation results of the application and the calculation results of the complex mirror image method and the 141-point linear filtering method is within 0.07, the error is very small, so the apparent resistivity calculation method of the horizontal layered earth of the application can effectively improve the calculation efficiency, and also has high calculation accuracy.

[0126] In the preferred embodiment of the present application, the integral variable set is obtained by the exponential discrete sampling method, the characteristic equation of the horizontally layered earth is constructed according to the integral variable set and the horizontal layered earth parameters of the target large area region obtained, the upper limit of the horizontally layered earth integral is obtained according to the first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation under each integral variable, the upper limit of the horizontally layered earth integral is set as the maximum value of the integral variable satisfying the preset criterion, the preset criterion is used to reflect the relationship between the local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set and the allowable error, the allowable error is determined according to the fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with the integral variable, the unbounded domain apparent resistivity relationship of the horizontally layered earth is constructed according to the pole distance of the Wenner quadrupole method, the integral range of the unbounded domain apparent resistivity relationship is limited based on the upper limit of the horizontally layered earth integral, the bounded domain apparent resistivity relationship is obtained, and the iterative solution of the bounded domain apparent resistivity relationship is performed based on the 151-point Gauss-Legendre integral algorithm and Euler transformation to obtain the apparent resistivity calculation result. The apparent resistivity calculation method of the horizontally layered earth disclosed in the present application determines the integral variable by using the equal exponential sampling method, can obtain the sharp change of the first layered uplink electromagnetic wave amplitude kernel function at a low frequency band, can also obtain the asymptotic value of the first layered uplink electromagnetic wave amplitude kernel function at a high frequency band, and thus can obtain the dynamic change waveform of the entire first layered uplink electromagnetic wave amplitude kernel function, thereby ensuring the accuracy of the calculation. At the same time, the integral variable exponential discretization step is determined, avoiding the influence of soil parameter changes. The upper limit of the horizontally layered earth integral is determined by the soil parameters, reflecting the difference of the soil parameters, and there is no need to use a large number of cases to solve the optimal solution of the filtering coefficient and the exponential sampling coefficient, thereby significantly reducing the complexity of the calculation. The Euler transformation of the apparent resistivity is realized by using the 151-point Gauss-Legendre numerical integral algorithm combined with the Euler transformation coefficient, which can effectively improve the calculation efficiency of the apparent resistivity of the horizontally layered earth and reduce the complexity of the calculation.

[0127] Correspondingly, as Figure 3 shown is an apparent resistivity calculation system of a horizontally layered earth, based on an apparent resistivity calculation method of a horizontally layered earth, the embodiment of the present application further provides an apparent resistivity calculation system of a horizontally layered earth, realizes the apparent resistivity calculation method of the horizontally layered earth disclosed in the embodiment of the present application, and includes a characteristic equation construction unit 1, a horizontally layered earth integral upper limit calculation unit 2, an apparent resistivity relationship construction unit 3, an apparent resistivity relationship conversion unit 4 and an apparent resistivity calculation unit 5.

[0128] The characteristic equation construction unit 1 is configured to obtain an integral variable set by an exponential discrete sampling method, and construct a characteristic equation of horizontally layered ground according to the integral variable set and horizontal layered ground parameters of a target large area region obtained;

[0129] The horizontally layered ground integral upper limit calculation unit 2 is configured to obtain a horizontally layered ground integral upper limit according to a first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation under each integral variable, the horizontally layered ground integral upper limit being set as a maximum value of the integral variable satisfying a preset criterion, the preset criterion being configured to reflect a relationship between a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set and an allowable error, the allowable error being set according to a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with the integral variable;

[0130] The apparent resistivity relationship construction unit 3 is configured to construct an unbounded domain apparent resistivity relationship of the horizontally layered ground according to a pole distance of the Wenner quadrupole method obtained;

[0131] The apparent resistivity relationship conversion unit 4 is configured to limit an integral range of the unbounded domain apparent resistivity relationship based on the horizontally layered ground integral upper limit, and obtain a bounded domain apparent resistivity relationship;

[0132] The apparent resistivity calculation unit 5 is configured to iteratively solve the bounded domain apparent resistivity relationship based on a 151-point Gauss-Legendre integral algorithm and Euler transformation, and obtain an apparent resistivity calculation result.

[0133] Further, the characteristic equation construction unit 1 includes a reflection coefficient calculation module, an interlayer relationship construction module, and a characteristic equation construction module;

[0134] The reflection coefficient calculation module is configured to obtain horizontal layered ground parameters of a target large area region, the horizontal layered ground parameters including soil resistivity of each layer and z-axis coordinates of layer boundaries of each layer, and construct a reflection coefficient of each layer boundary according to the soil resistivity;

[0135] The interlayer relationship construction module is configured to determine an interlayer relationship of layer boundaries of the horizontally layered ground according to the reflection coefficient, the integral variable set, and the z-axis coordinates of the layer boundaries based on electromagnetic characteristics of the horizontally layered ground;

[0136] The characteristic equation construction module is configured to construct a characteristic equation of the horizontally layered earth based on electromagnetic field boundary conditions, the characteristic equation being composed of a characteristic equation matrix, a characteristic equation vector and a to-be-solved eigenvector, elements of the characteristic equation matrix representing the interlayer relationship, elements of the characteristic equation vector representing known electromagnetic field boundary excitation conditions, and the to-be-solved eigenvector including the first layered uplink electromagnetic wave amplitude kernel function.

[0137] Further, the horizontally layered earth integral upper limit calculation unit 2 includes a first calculation module, an allowable error construction module, a local absolute error construction module and a second calculation module.

[0138] The first calculation module is configured to iteratively solve the characteristic equation according to the integral variable set to obtain the to-be-solved eigenvector corresponding to each integral variable in the integral variable set, and obtain the first layered uplink electromagnetic wave amplitude kernel function in each to-be-solved eigenvector.

[0139] The allowable error construction module is configured to obtain a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function according to each first layered uplink electromagnetic wave amplitude kernel function, and construct an allowable error according to the fluctuation range.

[0140] The local absolute error construction module is configured to obtain a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set according to each first layered uplink electromagnetic wave amplitude kernel function.

[0141] The second calculation module is configured to construct a preset criterion according to the local absolute error and the allowable error, and solve the preset criterion to obtain a maximum value of the integral variable that satisfies the preset criterion, and take the maximum value as the horizontally layered earth integral upper limit.

[0142] Further, the apparent resistivity calculation unit 3 includes a third calculation module, an Euler transformation module and a fourth calculation module.

[0143] The third calculation module is configured to solve the bounded domain apparent resistivity relationship by using a 151-point Gauss-Legendre integral algorithm to obtain an apparent resistivity Gauss calculation formula.

[0144] The Euler transformation module is configured to introduce an Euler transformation coefficient function into the apparent resistivity Gauss calculation formula to obtain an apparent resistivity Euler transformation calculation formula.

[0145] The fourth calculation module is configured to iteratively solve the apparent resistivity Euler transformation calculation formula to obtain an apparent resistivity measurement result.

[0146] Further, the apparent resistivity Euler transform calculation formula is expressed as:

[0147]

[0148] wherein, represents the pole distance, represents the pole distance number of the Wenner quadrupole method, represents the pole distance corresponding apparent resistivity, represents the soil resistivity of the first layer, represents the integral point index of the 151-point Gauss-Legendre integral, represents the position coefficient of the 151-point Gauss-Legendre integral, represents the upper limit of horizontal layered earth integral, represents the weight coefficient of the 151-point Gauss-Legendre integral, represents the coefficient function of Euler transform, represents the integral variable corresponding first-layer uplink electromagnetic wave amplitude kernel function, represents corresponding first-order zero Bessel function of the first kind, represents corresponding first-order zero Bessel function of the first kind

[0149] The specific limitations of the apparent resistivity measurement system of a horizontally layered earth can be referred to the above limitations of the apparent resistivity measurement method of a horizontally layered earth, which will not be repeated here. Those skilled in the art can realize that the various modules and steps described in combination with the embodiments disclosed in the present application can be realized in hardware, software or combination of both. Whether the functions are realized in hardware or software depends on the specific application and design constraints of the technical solution. The skilled person can use different methods to realize the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0150] In summary, the application provides a kind of apparent resistivity measurement method and system of horizontal layered earth, which solves the technical problem of how to meet the dual requirements of precision and efficiency of apparent resistivity measurement of horizontal layered earth, the method comprises: obtaining integral variable set by exponential discrete sampling method, and constructing the characteristic equation of horizontal layered earth according to the integral variable set and the horizontal layered earth parameters of the target large area region obtained;According to the first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation under each integral variable, the upper limit of horizontal layered earth integral is obtained, the upper limit of horizontal layered earth integral is set as the maximum value of integral variable that meets the preset criterion, the preset criterion is used to reflect the relationship between the local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set and the allowable error, and the allowable error is determined according to the fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function with integral variable;According to the pole distance of Wenner quadrupole method obtained, the unbounded domain apparent resistivity relationship of horizontal layered earth is constructed;Based on the upper limit of horizontal layered earth integral, the integral range of unbounded domain apparent resistivity relationship is limited to obtain bounded domain apparent resistivity relationship;Based on 151-point Gauss-Legendre integral algorithm and Euler transformation, the bounded domain apparent resistivity relationship is iteratively solved to obtain the apparent resistivity calculation result.The apparent resistivity measurement method of horizontal layered earth disclosed in the application uses equal exponential sampling method to determine the integral variable, which can obtain the sharp change of the first layered uplink electromagnetic wave amplitude kernel function in low frequency band, and can also obtain the asymptotic value of the first layered uplink electromagnetic wave amplitude kernel function in high frequency band, so as to obtain the dynamic change waveform of the entire first layered uplink electromagnetic wave amplitude kernel function, and ensure the accuracy of the calculation.At the same time, the integral variable exponential discretization step is determined, which avoids the influence of soil parameter change on it.The upper limit of horizontal layered earth integral is determined by soil parameters, which reflects the difference of soil parameters, and does not need to use a large number of cases to solve the optimal solution of filtering coefficient and exponential sampling coefficient, which significantly reduces the complexity of calculation.Using 151-point Gauss-Legendre numerical integral algorithm combined with Euler transformation coefficient to realize the Euler transformation of apparent resistivity can effectively improve the calculation efficiency of apparent resistivity of horizontal layered earth and reduce the complexity of calculation.

[0151] The various embodiments in this specification are described in a progressive manner. For directly identical or similar parts of the embodiments, refer to each other. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions in the method embodiments. It should be noted that the technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification.

[0152] The embodiments described above are merely preferred embodiments of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the patent application. It should be noted that those skilled in the art can make various improvements and substitutions without departing from the technical principles of this application, and these improvements and substitutions should also be considered within the scope of protection of this application. Therefore, the scope of protection of this patent application should be determined by the scope of the claims.

Claims

1. A method of measuring apparent resistivity of a horizontally layered earth, characterized by, The method comprises: obtaining an integral variable set by an exponential discrete sampling method, and constructing a characteristic equation of horizontally layered earth according to the integral variable set and acquired horizontal layered earth parameters of a target large area; obtaining a horizontally layered earth integral upper limit according to a first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation under each integral variable, the horizontally layered earth integral upper limit being set as a maximum value of the integral variable satisfying a preset criterion, the preset criterion being used to reflect a relationship between a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to a highest integral variable in the integral variable set and an allowable error, the allowable error being set as a fluctuation range determined according to the first layered uplink electromagnetic wave amplitude kernel function varying with the integral variable; constructing a non-bounded domain apparent resistivity relationship of the horizontally layered earth according to an acquired pole distance of the Wenner quadrupole method; limiting an integral range of the non-bounded domain apparent resistivity relationship based on the horizontally layered earth integral upper limit to obtain a bounded domain apparent resistivity relationship; iteratively solving the bounded domain apparent resistivity relationship based on a 151-point Gauss-Legendre integral algorithm and Euler transformation to obtain an apparent resistivity calculation result.

2. The method of apparent resistivity sounding of horizontally layered earth as claimed in claim 1 wherein, The method comprises: acquiring horizontal layered earth parameters of a target large area, the horizontal layered earth parameters comprising soil resistivity of each layer and z-axis coordinates of layer boundaries of each layer, and constructing a reflection coefficient of each layer boundary according to the soil resistivity; determining an interlayer relationship of layer boundaries of the horizontally layered earth according to the reflection coefficient, the integral variable set and the z-axis coordinates of the layer boundaries based on electromagnetic characteristics of the horizontally layered earth; constructing a characteristic equation of the horizontally layered earth based on electromagnetic field boundary conditions, the characteristic equation being composed of a characteristic equation matrix, a characteristic equation vector and a to-be-solved characteristic vector, elements of the characteristic equation matrix representing the interlayer relationship, elements of the characteristic equation vector representing known electromagnetic field boundary excitation conditions, and the to-be-solved characteristic vector comprising the first layered uplink electromagnetic wave amplitude kernel function.

3. The method of apparent resistivity sounding of horizontally layered earth as claimed in claim 2 wherein, The method comprises: iteratively solving the characteristic equation according to the integral variable set to obtain the to-be-solved characteristic vector corresponding to each integral variable in the integral variable set, and acquiring the first layered uplink electromagnetic wave amplitude kernel function in each to-be-solved characteristic vector; obtaining a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function according to each first layered uplink electromagnetic wave amplitude kernel function, and constructing an allowable error according to the fluctuation range; obtaining a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to a highest integral variable in the integral variable set according to each first layered uplink electromagnetic wave amplitude kernel function. According to the local absolute error and the allowed error, a preset criterion is constructed, and the preset criterion is solved to obtain a maximum value of the integral variable satisfying the preset criterion, and the maximum value is taken as a horizontal layered geodetic integral upper limit.

4. The method of apparent resistivity sounding of horizontally layered earth as claimed in claim 1 wherein, The 151-point Gauss-Legendre integral algorithm and Euler transformation are used to iteratively solve the bounded domain apparent resistivity relationship to obtain an apparent resistivity calculation result, including: The 151-point Gauss-Legendre integral algorithm is used to solve the bounded domain apparent resistivity relationship to obtain an apparent resistivity Gauss calculation formula; An Euler transformation coefficient function is introduced into the apparent resistivity Gauss calculation formula to obtain an apparent resistivity Euler transformation calculation formula; The apparent resistivity Euler transformation calculation formula is iteratively solved to obtain an apparent resistivity calculation result.

5. The method of apparent resistivity sounding of horizontally layered earth as claimed in claim 4 wherein, The apparent resistivity Euler transformation calculation formula is expressed as: wherein, denotes the pole distance, denotes the pole distance number of the Weyner quadrupole method, denotes the pole distance corresponding apparent resistivity, denotes the soil resistivity of the first layer, denotes the integration point index of the 151-point Gauss-Legendre integration, denotes the position coefficient of the 151-point Gauss-Legendre integration, denotes the upper limit of the horizontal layered earth integration, denotes the weight coefficient of the 151-point Gauss-Legendre integration, denotes the coefficient function of the Euler transformation, denotes the integration variable is corresponding first-layer uplink electromagnetic wave amplitude kernel function, denotes corresponding first type zero order Bessel function, denotes corresponding first type zero order Bessel function.

6. A system for measuring apparent resistivity of horizontally layered earth for implementing the method for measuring apparent resistivity of horizontally layered earth according to any one of claims 1 to 5, characterized in that, The system comprises a characteristic equation construction unit, a horizontal layered geodetic integral upper limit calculation unit, an apparent resistivity relationship construction unit, an apparent resistivity relationship conversion unit, and an apparent resistivity calculation unit; The characteristic equation construction unit is configured to obtain an integral variable set by an exponential discrete sampling method, and construct a characteristic equation of horizontal layered earth according to the integral variable set and acquired horizontal layered geodetic parameters of a target large area; The horizontal layered geodetic integral upper limit calculation unit is configured to obtain a horizontal layered geodetic integral upper limit according to a first layered uplink electromagnetic wave amplitude kernel function corresponding to the characteristic equation of each integral variable, wherein the horizontal layered geodetic integral upper limit is set as a maximum value of the integral variable satisfying a preset criterion, and the preset criterion is used to reflect a relationship between a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to a highest integral variable in the integral variable set and an allowed error determined according to a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function varying with the integral variable; The apparent resistivity relationship construction unit is configured to construct a unbounded domain apparent resistivity relationship of the horizontal layered earth according to an acquired pole distance of the Wenner quadrupole method; The apparent resistivity relationship conversion unit is configured to limit an integral range of the unbounded domain apparent resistivity relationship based on the horizontal layered geodetic integral upper limit to obtain a bounded domain apparent resistivity relationship. The apparent resistivity calculation unit is configured to iteratively solve the bounded domain apparent resistivity relationship based on the 151-point Gauss-Legendre integral algorithm and Euler transformation to obtain an apparent resistivity calculation result.

7. The system for measuring apparent resistivity of horizontally layered earth according to claim 6, wherein, The characteristic equation construction unit comprises a reflection coefficient calculation module, a layer relationship construction module, and a characteristic equation construction module. The reflection coefficient calculation module is configured to acquire horizontal layered geodetic parameters of a target large area, wherein the horizontal layered geodetic parameters comprise a soil resistivity of each layer and a layer boundary z-axis coordinate of each layer, and construct a reflection coefficient of each layer boundary according to the soil resistivity. The interlayer relationship construction module is configured to determine an interlayer relationship of a layered boundary of the horizontally layered earth based on electromagnetic characteristics of the horizontally layered earth, according to the reflection coefficient, the integral variable set, and the z-axis coordinate of the layered boundary. The characteristic equation construction module is configured to construct a characteristic equation of the horizontally layered earth based on electromagnetic field boundary conditions, the characteristic equation being composed of a characteristic equation matrix, a characteristic equation vector, and a to-be-solved characteristic vector, elements of the characteristic equation matrix representing the interlayer relationship, elements of the characteristic equation vector representing known electromagnetic field boundary excitation conditions, and the to-be-solved characteristic vector including the first layered uplink electromagnetic wave amplitude kernel function.

8. The system for measuring apparent resistivity of horizontally layered earth according to claim 7, wherein, The horizontally layered earth integral upper limit calculation unit includes a first calculation module, an allowable error construction module, a local absolute error construction module, and a second calculation module. The first calculation module is configured to iteratively solve the characteristic equation according to the integral variable set, to obtain the to-be-solved characteristic vector corresponding to each integral variable in the integral variable set, and to obtain the first layered uplink electromagnetic wave amplitude kernel function in each to-be-solved characteristic vector. The allowable error construction module is configured to obtain a fluctuation range of the first layered uplink electromagnetic wave amplitude kernel function according to each first layered uplink electromagnetic wave amplitude kernel function, and to construct an allowable error according to the fluctuation range. The local absolute error construction module is configured to obtain a local absolute error of the first layered uplink electromagnetic wave amplitude kernel function corresponding to the highest integral variable in the integral variable set according to each first layered uplink electromagnetic wave amplitude kernel function. The second calculation module is configured to construct a preset criterion according to the local absolute error and the allowable error, to solve the preset criterion, to obtain a maximum value of the integral variable that satisfies the preset criterion, and to take the maximum value as the horizontally layered earth integral upper limit.

9. The system for measuring apparent resistivity of horizontally layered earth according to claim 6, wherein, The apparent resistivity calculation unit includes a third calculation module, an Euler transformation module, and a fourth calculation module. The third calculation module is configured to solve the bounded domain apparent resistivity relationship using a 151-point Gauss-Legendre integral algorithm to obtain an apparent resistivity Gauss calculation formula. The Euler transformation module is configured to introduce an Euler transformation coefficient function into the apparent resistivity Gauss calculation formula to obtain an apparent resistivity Euler transformation calculation formula. The fourth calculation module is configured to iteratively solve the apparent resistivity Euler transformation calculation formula to obtain an apparent resistivity measurement result.

10. The system for measuring apparent resistivity of horizontally layered earth according to claim 9, wherein, The apparent resistivity Euler transformation calculation formula is represented as: wherein, denotes the pole distance, denotes the pole distance number of the Weyner quadrupole method, denotes the pole distance the corresponding apparent resistivity, denotes the soil resistivity of the first layer, denotes the integration point index of the 151-point Gauss-Legendre integration, denotes the position coefficient of the 151-point Gauss-Legendre integration, denotes the upper limit of the horizontal layered earth integration, denotes the weight coefficient of the 151-point Gauss-Legendre integration, denotes the coefficient function of the Euler transformation, denotes the integration variable is the corresponding first-layer uplink electromagnetic wave amplitude kernel function, denotes the corresponding first-order zero Bessel function of the first kind, denotes the corresponding first-order zero Bessel function of the first kind.

Images