A horizontal layered earth parameter inversion method, system, device and medium
By constructing the apparent resistivity prediction function and using odd polynomials, even polynomials and trust region algorithm, the problem of low precision in the inversion of horizontally layered geodetic parameters is solved, and the efficiency and accuracy of the inversion results are achieved.
Patent Information
- Application Number
- CN202511073817.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-08-01
- Publication Date
- 2025-10-10
- Estimated Expiration
- 2045-08-01
AI Technical Summary
In the existing technology, the accuracy of horizontal layered geodetic parameter inversion is low, which is difficult to meet the needs of high-precision engineering applications and geological analysis. This is mainly due to the problems of computational error accumulation and oscillation in the iterative process of deterministic optimization algorithms.
The apparent resistivity prediction function is constructed using kernel functions and first-kind zero-order Bessel functions. Odd and even polynomials are combined, and the partial derivative relationship is calculated using the permutation and combination algorithm and the complex mirror method. The trust region algorithm is used for inversion, which simplifies the response calculation of multi-layer soil media and improves the calculation accuracy and stability of the partial derivatives.
The accuracy and reliability of horizontal layered geodetic parameter inversion are improved, the efficiency and precision of the inversion results are ensured, and the oscillation and divergence problems in numerical calculations are avoided.
Smart Images

Figure CN120577884B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of geodetic parameter inversion, and in particular to a horizontal layered geodetic parameter inversion method, system, equipment and medium. Background Art
[0002] Accurately acquiring layered earth and soil model parameters is crucial for engineering design and geological analysis in fields such as power system grounding design, geological exploration, and electromagnetic environment monitoring. Horizontal layered earth parameter inversion is the core technical means for accurately acquiring these parameters. This inversion involves deducing a layered earth and soil model from apparent resistivity measurements taken in the field.
[0003] At present, deterministic optimization algorithms in horizontally layered geodetic parameter inversion are widely used in inversion calculations due to their clear iteration rules and fast local convergence speed. However, the convergence efficiency of deterministic optimization algorithms is highly dependent on the precise partial derivative data of the inversion objective function. The partial derivative data guides the parameter update direction, thereby accelerating the inversion process.
[0004] However, the theoretical derivation of the partial derivatives of the inversion objective function with respect to the soil resistivity of the layered earth and the inversion objective function with respect to the layer thickness is extremely complex. The existing technology usually adopts the numerical difference method for approximate calculation, but this method has significant defects. The selection of the difference step size directly affects the approximate accuracy of the partial derivative. Too large a step size can easily lead to data distortion, and too small a step size will introduce serious numerical noise, aggravating the calculation error. Rounding errors and truncation errors continue to accumulate during the calculation process, reducing the stability of the calculation method, causing the deterministic optimization algorithm to oscillate or even diverge during the iteration process, seriously restricting the calculation accuracy of the deterministic optimization algorithm, resulting in low accuracy of the inversion results of horizontally layered earth parameters, which is difficult to meet the needs of high-precision engineering applications and geological analysis.
[0005] It can be seen that how to improve the accuracy of horizontal layered geodetic parameter inversion has become a technical problem that needs to be solved urgently by those skilled in the art. Summary of the Invention
[0006] The present invention provides a horizontal layered geodetic parameter inversion method, system, device and medium to solve the technical problem of how to improve the accuracy of horizontal layered geodetic parameter inversion, thereby achieving the effect of improving the accuracy and reliability of the horizontal layered geodetic parameter inversion results.
[0007] In a first aspect, the present invention provides a method for inverting horizontally layered earth parameters, the method comprising:
[0008] Obtaining apparent resistivity measurement values and corresponding pole distances of a target earth region, and constructing an inversion objective function based on an apparent resistivity prediction function corresponding to the pole distance constructed using a kernel function and a first-kind zero-order Bessel function and the apparent resistivity measurement values;
[0009] The kernel function is represented by odd polynomials and even polynomials, and a permutation and combination algorithm is used to obtain a plurality of ordered sets. The odd polynomials and the even polynomials are represented by the plurality of ordered sets and then partial derivatives are calculated to obtain a first partial derivative relationship between the kernel function and soil resistivity and a second partial derivative relationship between the kernel function and layer thickness, wherein the ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums;
[0010] According to the first partial derivative relationship and the second partial derivative relationship, the first partial derivative calculation formula and the second partial derivative calculation formula are respectively subjected to complex mirror fitting using a complex mirror method to obtain a third partial derivative relationship and a fourth partial derivative relationship, wherein the first partial derivative calculation formula and the third partial derivative relationship are both used for calculating the partial derivative of the inversion objective function with respect to the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for calculating the partial derivative of the inversion objective function with respect to the layer thickness;
[0011] According to the third partial derivative relationship and the fourth partial derivative relationship, a trust region algorithm is used to perform horizontal layered earth parameter inversion to obtain soil resistivity values and layer thickness values of each layer in the target earth region.
[0012] Preferably, the permutation and combination algorithm is used to obtain a plurality of ordered sets, and the odd polynomials and the even polynomials are represented by the plurality of ordered sets and then partial derivative calculations are performed to obtain a first partial derivative relationship of the kernel function to soil resistivity and a second partial derivative relationship of the kernel function to layer thickness, including:
[0013] Constructing a first set according to the soil layer number, setting the number of iterations, and using a permutation and combination algorithm to arbitrarily select the number of elements from the first set and arrange them in ascending order to obtain a second set;
[0014] Using the elements of the second set as numbers, and obtaining the soil resistivity change rate of adjacent layers corresponding to each number according to a preset first calculation rule, permuting and combining the numbers of the soil resistivity change rates of adjacent layers in ascending order to obtain the soil resistivity change rate set of adjacent layers, obtaining the layer thickness sum corresponding to each number according to a preset second calculation rule, permuting and combining the numbers of the layer thickness sums in ascending order to obtain the layer thickness sum set, obtaining the segmented value corresponding to each number according to a preset segmented value assignment rule, and permuting and combining the numbers of the segmented value values in ascending order to obtain a third set;
[0015] Based on the set of soil resistivity change rates of adjacent layers and the set of layer thickness sums, the odd polynomial is represented according to the properties of the odd polynomial and then iteratively calculated to obtain an odd polynomial expression; and the even polynomial is represented according to the properties of the even polynomial and then iteratively calculated to obtain an even polynomial expression;
[0016] Obtaining a third partial derivative calculation formula of the kernel function with respect to the soil resistivity, and updating the third partial derivative calculation formula according to the odd polynomial expression and the even polynomial expression to obtain a first partial derivative relationship formula of the kernel function with respect to the soil resistivity;
[0017] The fourth partial derivative calculation formula of the kernel function for the layer thickness is updated according to the odd polynomial expression, the even polynomial expression and the third set to obtain the second partial derivative relationship of the kernel function for the layer thickness.
[0018] Preferably, the odd polynomial expression is:
[0019]
[0020] The even polynomial expression is:
[0021]
[0022] in, represents the integration variable, represents the number of iterations, represents the total number of soil stratification layers, represents the elements of the second set, represents the first element of the second set, represents the second element of the second set, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the numbered layer thickness and, represents the set of soil resistivity change rates of adjacent layers, represents the layer thickness and set, Indicates from Choose from natural numbers The number of combinations of natural numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the numbered layer thickness and, Represents the elements of the second set is the thickness and of the numbered layers.
[0023] Preferably, the apparent resistivity prediction function is expressed as:
[0024]
[0025] in, Indicates the The measured pole distance, Indicates pole distance The corresponding predicted apparent resistivity value is, Indicates the Layered soil resistivity, Indicates the Hierarchical Hankel integral kernel function, express The corresponding zero-order Bessel function of the first kind is, express The corresponding zero-order Bessel function of the first kind is, represents the integration variable.
[0026] Preferably, constructing an inversion objective function based on the apparent resistivity prediction function corresponding to the pole distance constructed using a kernel function and a first-kind zero-order Bessel function and the apparent resistivity measurement value comprises:
[0027] Obtaining a measurement error according to the apparent resistivity prediction function and the apparent resistivity measurement value;
[0028] An inversion objective function is constructed with the goal of minimizing the ratio of the measurement error to the apparent resistivity measurement value.
[0029] Preferably, the first partial derivative calculation formula and the second partial derivative calculation formula are respectively subjected to complex mirror fitting by the complex mirror method according to the first partial derivative relationship and the second partial derivative relationship, to obtain a third partial derivative relationship and a fourth partial derivative relationship, respectively, including:
[0030] Obtaining a first kernel function correlation formula in a fifth partial derivative calculation formula of the apparent resistivity prediction function with respect to the soil resistivity, and a second kernel function correlation formula in a sixth partial derivative calculation formula of the apparent resistivity prediction function with respect to the layer thickness;
[0031] According to the first partial derivative relationship and the second partial derivative relationship, a complex mirror method is used to perform complex mirror fitting on the first kernel function correlation relationship and the second kernel function correlation relationship, respectively, to obtain a first complex mirror fitting result and a second complex mirror fitting result;
[0032] Obtaining a first partial derivative calculation formula of the inversion objective function with respect to the soil resistivity, and sequentially updating the fifth partial derivative calculation formula and the first partial derivative calculation formula according to the first complex mirror fitting result to obtain a third partial derivative relationship formula of the inversion objective function with respect to the soil resistivity;
[0033] Obtain a second partial derivative calculation formula of the inversion objective function for the layer thickness, and update the sixth partial derivative calculation formula and the second partial derivative calculation formula in sequence according to the second complex mirror fitting result to obtain a fourth partial derivative relationship of the inversion objective function for the layer thickness.
[0034] Preferably, the inversion of horizontally layered earth parameters using a trust region algorithm based on the third partial derivative relationship and the fourth partial derivative relationship to obtain soil resistivity values and layer thickness values of each layer in the target earth region includes:
[0035] Determine the inversion starting point, and obtain the soil resistivity partial derivative of the inversion objective function for each layer according to the third partial derivative relationship, and obtain the layer thickness partial derivative of the inversion objective function for each layer according to the fourth partial derivative relationship;
[0036] Obtaining an initial gradient vector of an inversion objective function according to the partial derivative of the soil resistivity and the partial derivative of the layer thickness corresponding to each layer, and setting an upper limit of a trust region radius and an initial trust region radius according to the initial gradient vector of the objective function;
[0037] During the iterative inversion process, a current gradient vector of the inversion objective function is calculated, and it is determined whether the current gradient vector of the inversion objective function satisfies a preset convergence judgment condition. If so, an inversion result is output, wherein the inversion result includes soil resistivity values and soil thickness values of each layer of the target earth region;
[0038] If not, solving the optimal step length within the current trust region radius limit according to the third partial derivative relationship and the fourth partial derivative relationship, and calculating the ratio of the actual descent amount of the inversion objective function to the predicted descent amount of the inversion objective function according to the optimal step length, and iteratively correcting the current trust region radius according to the ratio;
[0039] Determining whether the ratio is greater than a preset lower limit proportional coefficient, if so, updating the current inversion point and updating the feasible search direction matrix based on the updated current inversion point, if not, maintaining the current inversion point;
[0040] The number of iterations is updated to return to the iterative inversion process to continue iterative operations until the inversion result is output.
[0041] In a second aspect, the present invention further provides a horizontally layered geodetic parameter inversion system for implementing the above-mentioned horizontally layered geodetic parameter inversion method, the system comprising: an inversion objective function construction module, a kernel function partial derivative deduction module, an inversion objective function partial derivative deduction module, and a horizontally layered geodetic parameter inversion module;
[0042] The inversion objective function construction module is used to obtain the apparent resistivity measurement value and the corresponding pole distance of the target earth area, and construct the inversion objective function according to the apparent resistivity prediction function corresponding to the pole distance constructed using the kernel function and the first-kind zero-order Bessel function and the apparent resistivity measurement value;
[0043] The kernel function partial derivative deduction module is used to represent the kernel function using odd polynomials and even polynomials, and to obtain a plurality of ordered sets using a permutation and combination algorithm. The odd polynomials and the even polynomials are represented using the plurality of ordered sets and then partial derivative calculations are performed to obtain a first partial derivative relationship between the kernel function and soil resistivity and a second partial derivative relationship between the kernel function and layer thickness. The ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums.
[0044] The inversion objective function partial derivative deduction module is used to perform complex mirror fitting on the first partial derivative calculation formula and the second partial derivative calculation formula respectively according to the first partial derivative relationship and the second partial derivative relationship using the complex mirror method, and obtain a third partial derivative relationship and a fourth partial derivative relationship respectively. The first partial derivative calculation formula and the third partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the layer thickness;
[0045] The horizontal layered earth parameter inversion module is used to perform horizontal layered earth parameter inversion based on the third partial derivative relationship and the fourth partial derivative relationship using a trust region algorithm to obtain the soil resistivity value and layer thickness value of each layer in the target earth area.
[0046] In a third aspect, the present invention also provides a computer device, which includes a memory, a processor and a transceiver, which are connected via a bus; the memory is used to store a set of computer program instructions and data, and transmit the stored data to the processor, and the processor executes the computer program instructions stored in the memory to execute the above-mentioned horizontal layered earth parameter inversion method.
[0047] In a fourth aspect, the present invention further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed, the above-mentioned horizontal layered earth parameter inversion method is implemented.
[0048] This application provides a method, system, device, and medium for inverting horizontally layered geodetic parameters. Compared with the prior art, the embodiments of this application have the following beneficial effects:
[0049] The horizontal layered geodetic parameter inversion method disclosed in this application uses the Hankel integral kernel function to describe the axisymmetric field problem of radial diffusion of current in the stratum in the four-probe method, converting the complex spatial domain problem into algebraic operations in the wavenumber domain, greatly simplifying the response calculation of multi-layered soil media. The smooth mathematical properties of the first-kind zero-order Bessel function avoid oscillation or divergence problems in numerical calculations, maintain the stability of the apparent resistivity prediction function, and improve the accuracy and reliability of the inversion results. The introduction of odd and even polynomials, combined with the permutation and combination algorithm, effectively improves the calculation speed and accuracy of the partial derivative of the inversion objective function with respect to soil resistivity and the partial derivative of the inversion objective function with respect to layer thickness, ensuring the high efficiency of horizontal layered geodetic parameter inversion in the target geodetic region and the accuracy of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS
[0050] Figure 1 This is a schematic diagram of the steps of a horizontal layered earth parameter inversion method provided by a preferred embodiment of the present invention;
[0051] Figure 2 This is a schematic structural diagram of a horizontally layered earth parameter inversion system provided by a preferred embodiment of the present invention;
[0052] Figure 3 is an internal structural diagram of a computer device according to an embodiment of the present invention;
[0053] Reference numerals:
[0054] 1- Inversion objective function construction module, 2- Kernel function partial derivative deduction module, 3- Inversion objective function partial derivative deduction module, 4- Horizontal layered earth parameter inversion module. DETAILED DESCRIPTION
[0055] The following is a detailed explanation of the embodiments of the present invention in conjunction with the accompanying drawings. The embodiments are provided for illustrative purposes only and cannot be understood as limitations on the present invention. The accompanying drawings are for reference and illustration purposes only and do not constitute a limitation on the scope of patent protection of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. In the description of the present invention, the terms "first", "second", "third", etc. are only used for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined as "first", "second", "third", etc. may explicitly or implicitly include one or more of the features. In the description of the present invention, unless otherwise specified, the meaning of "multiple" is two or more.
[0056] In the description of the present invention, it should be noted that, unless otherwise expressly specified and limited, the terms "installed", "connected" and "connected" should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or an indirect connection through an intermediate medium, or it can be a communication between the two components. 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 component referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as a limitation on the present invention. The term "and / or" used herein includes any and all combinations of one or more related listed items. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.
[0057] In describing the present invention, it should be noted that, unless otherwise defined, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art. The terms used in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. Those skilled in the art will understand the specific meanings of the above terms in the present invention in specific circumstances.
[0058] See also Figure 1 FIG2 is a schematic diagram of a method for inverting horizontally layered earth parameters. In an embodiment of the present invention, a method for inverting horizontally layered earth parameters is provided, and the method includes:
[0059] S1. Obtain the apparent resistivity measurement value and the corresponding pole distance of the target earth area, and construct an inversion target function based on the apparent resistivity prediction function corresponding to the pole distance constructed using a kernel function and a first-order zero-order Bessel function and the apparent resistivity measurement value; in a preferred embodiment of the present application, a four-probe method is used to obtain the apparent resistivity measurement value of the target earth area. The four-probe method includes four electrodes, two power supply electrodes and two measuring electrodes. Current is applied through the power supply electrode and the potential difference between the power supply electrode and the measuring electrode is collected through the measuring electrode to calculate the apparent resistivity of the target earth area. The distance between the power supply electrode and the measuring electrode is the pole distance. Set the horizontal M-layer earth parameters used for the inversion of the horizontal layered earth parameters of the target earth area, including the soil resistivity, layer thickness and layer boundary z-axis coordinate corresponding to each layer. During the four-probe measurement process, the input items include N groups of apparent resistivity measurement data, which are expressed as follows:
[0060]
[0061] in, Indicates the number of measurements, represents the total number of measurements, Indicates the The apparent resistivity value of the measurement, Indicates the The larger the pole distance, the deeper the detection depth. Therefore, it is necessary to ensure that the maximum pole distance covers the maximum layer depth of the horizontal layered geodetic parameter inversion model, and the number of measured data points must not be less than the number of inversion layers to avoid underdetermination.
[0062] The horizontal M-layer geodetic parameters used in the horizontal layered geodetic parameter inversion model include the soil resistivity, layer thickness, and layer boundary z-axis coordinate corresponding to each layer, which can be specifically expressed as:
[0063]
[0064] in, Indicates the layer number, Indicates the Layered soil resistivity, Indicates the Layer thickness of the layers, Indicates the The z-axis coordinate of the layered bottom surface, Indicates the The z-coordinate of the bottom layer.
[0065] In a preferred embodiment of the present application, a kernel function and a first-kind zero-order Bessel function are used to construct an apparent resistivity prediction function corresponding to the pole distance. The Hankel integral kernel function is used as the kernel function to transform the partial differential equation in the spatial domain into an algebraic equation in the wavenumber domain through an integral transformation, thereby simplifying the solution process. The first-kind zero-order Bessel function is a special solution of the Bessel equation when the order is equal to 0. The apparent resistivity prediction function is expressed as:
[0066]
[0067] in, Indicates the The measured pole distance, Indicates pole distance The corresponding predicted apparent resistivity value is, represents the soil resistivity of the first layer, represents the Hankel integral kernel function of the first layer, express The corresponding zero-order Bessel function of the first kind is, express The corresponding zero-order Bessel function of the first kind is, represents the integration variable.
[0068] In a preferred embodiment of the present application, a Hankel integral kernel function is used to describe the axisymmetric field problem of radial current diffusion in the stratum using the four-probe method. In a horizontally layered earth model, the resistivity distribution of the underground medium is uniform along the horizontal direction, meeting the cylindrical symmetry condition. The Hankel integral kernel function serves as a transformation kernel to relate the physical quantities of radial distance and vertical depth, transforming complex spatial domain problems into algebraic operations in the wavenumber domain, greatly simplifying the response calculation of multi-layered soil media. The smooth mathematical properties of the first-kind zero-order Bessel function avoid oscillation or divergence problems in numerical calculations, especially when processing large-pole-spacing data. It can maintain the stability of the apparent resistivity prediction function and improve the accuracy and reliability of the inversion results.
[0069] In a preferred embodiment of the present application, the measurement error is obtained based on the apparent resistivity prediction function and the apparent resistivity measurement value. With the goal of minimizing the ratio of the measurement error to the apparent resistivity measurement value, an inversion objective function is constructed. The inversion objective function is expressed as:
[0070]
[0071] in, represents the inversion objective function, Indicates the The apparent resistivity measurement value of the measurement.
[0072] S2. Use odd polynomials and even polynomials to represent the kernel function, use a permutation and combination algorithm to obtain several ordered sets, and use several of the ordered sets to represent the odd polynomials and the even polynomials and then perform partial derivative calculations, correspondingly obtaining the first partial derivative relationship of the kernel function to soil resistivity and the second partial derivative relationship of the kernel function to layer thickness, the ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums; in a preferred embodiment of the present application, the partial derivatives of the inversion objective function to the soil resistivity or layer thickness of each layer in the target earth area are converted into partial derivatives of the kernel function to the soil resistivity or layer thickness, and the permutation and combination algorithm is further used to derive a general expression of the kernel function in the apparent resistivity prediction function, and based on the general expression, the first partial derivative relationship of the kernel function to soil resistivity and the second partial derivative relationship of the kernel function to layer thickness are obtained.
[0073] According to the expression of the inversion objective function, the calculation formula for the first partial derivative of the inversion objective function with respect to soil resistivity is obtained as follows:
[0074]
[0075] The calculation formula for the second partial derivative of the inversion objective function with respect to the layer thickness is:
[0076]
[0077] According to the expression of the apparent resistivity prediction function, the calculation formula of the fifth partial derivative of the apparent resistivity prediction function with respect to soil resistivity is obtained, which is expressed as:
[0078]
[0079] The calculation formula of the sixth partial derivative of the apparent resistivity prediction function with respect to the layer thickness is expressed as:
[0080]
[0081] In the preferred embodiment of the present application, odd polynomials and even polynomials are introduced, and odd polynomials and even polynomials are used to calculate the first The layered Hankel integral kernel function is described as follows:
[0082]
[0083] in, represents an odd polynomial, Represents an even polynomial.
[0084] From the partial derivative formula of the apparent resistivity prediction function with respect to soil resistivity or the partial derivative formula of the apparent resistivity prediction function with respect to layer thickness, it can be seen that once the first partial derivative relationship of the kernel function with respect to soil resistivity and the second partial derivative relationship of the kernel function with respect to layer thickness are obtained, the third partial derivative relationship of the inversion objective function with respect to soil resistivity and the fourth partial derivative relationship of the inversion objective function with respect to layer thickness can be obtained using the complex mirror method. In this application, a permutation and combination algorithm is used to obtain several ordered sets. After representing odd and even polynomials based on these ordered sets and then performing partial derivative calculations, the first partial derivative relationship of the kernel function with respect to soil resistivity and the second partial derivative relationship of the kernel function with respect to layer thickness can be obtained.
[0085] Specifically, for the horizontal M layers of earth, the first set is constructed according to the number of soil layers, which is expressed as , set the number of iterations K, the initial value of the number of iterations is 0, the maximum value is M-1, when the number of iterations K is 0, =1, =0.
[0086] make , update the number of iterations, use the permutation and combination algorithm to select K elements from the first set in sequence and arrange them in ascending order to obtain the second set, which is expressed as , because the number of iterations K is less than or equal to M-1, the result of the permutation and combination is not unique. For example, when M = 5 and K = 2, the first set is , take 2 elements from the first set and arrange them in descending order, the second set is 、 、 、 、 、 , the elements in the second set correspond to the numbers of the horizontal layers.
[0087] Furthermore, the elements of the second set are used as numbers, and the adjacent layer soil resistivity change rate corresponding to each number is obtained according to a preset first calculation rule. The numbers of the adjacent layer soil resistivity change rates are arranged and combined in ascending order to obtain an adjacent layer soil resistivity change rate set. The first calculation rule is:
[0088]
[0089] in, Indicates the number The change rate of soil resistivity of adjacent layers, Indicates the Layered soil resistivity, Indicates the Layered soil resistivity.
[0090] The numbers of the calculated soil resistivity change rates of adjacent layers are arranged and combined in ascending order to obtain the set of soil resistivity change rates of adjacent layers, which is expressed as For M = 5, K = 2, the corresponding set of soil resistivity change rates of adjacent layers is 、 、 、 、 、 .
[0091] Furthermore, the layer thickness sum corresponding to each number is obtained according to a preset second calculation rule, and the layer thickness sums are arranged and combined to obtain a layer thickness sum set. The second calculation rule is:
[0092]
[0093] in, Indicates the number The layer thickness and Indicates the Layer thickness of the layer.
[0094] The calculated layer thickness and are arranged and combined to obtain the layer thickness and set, which is expressed as For M = 5, K = 2, the corresponding layer thickness and set are 、 、 、 、 、 .
[0095] Furthermore, the partial derivative algorithm is used to calculate the partial derivatives of the Hankel integral kernel function of the first layer represented by odd polynomials and even polynomials, and the third partial derivative calculation formula of the Hankel integral kernel function of the first layer to the soil resistivity and the fourth partial derivative calculation formula of the Hankel integral kernel function of the first layer to the layer thickness are obtained. The third partial derivative calculation formula is:
[0096]
[0097] The fourth partial derivative calculation formula is:
[0098]
[0099] When K is an even number, the operation formula for the even polynomial is:
[0100]
[0101] in, represents the number of iterations, represents the elements of the second set, represents the first element of the second set, represents the second element of the second set, represents the rate of change of soil resistivity of adjacent layers numbered by the elements of the second set, represents the sum of the thickness of the layers numbered by the elements of the second set, represents the set of soil resistivity change rates of adjacent layers, represents the layer thickness and set, Indicates from Choose from natural numbers The number of combinations of natural numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the numbered layer thickness and, Represents the elements of the second set is the thickness and of the numbered layers.
[0102] When K is an odd number, the operation formula for the odd polynomial is:
[0103]
[0104] After traversing the adjacent layer soil resistivity change rate set and layer thickness sum set, iterate the number of iterations K until , we get the odd polynomial expression and the even polynomial expression. The even polynomial expression is expressed as:
[0105]
[0106] The odd polynomial expression is expressed as:
[0107]
[0108] According to the first calculation rule, the partial derivative of the rate of change of soil resistivity of adjacent layers with respect to soil resistivity is calculated. Since the elements in the second set correspond to the numbers of the horizontal layers, the partial derivative of the rate of change of soil resistivity of adjacent layers with respect to soil resistivity can be expressed as:
[0109]
[0110]
[0111] After substituting the partial derivative of the soil resistivity change rate of adjacent layers with respect to the soil resistivity into the third partial derivative calculation formula, the updated third partial derivative calculation formula is obtained:
[0112]
[0113] Substituting the odd polynomial expression and the even polynomial expression into the updated third partial derivative calculation formula, the first partial derivative relationship of the kernel function to soil resistivity can be obtained.
[0114] For the second partial derivative of the Hankel integral kernel function of the first layer with respect to the layer thickness, it is necessary to obtain the segmented value corresponding to each number according to the preset segmented assignment rule, and to arrange and combine the segmented assignment numbers in ascending order to obtain the third set. The preset segmented assignment rule is expressed as:
[0115]
[0116] The third set is represented as For M = 5, K = 2, the corresponding layer thickness and set are 、 、 、 、 、 .
[0117] According to the third set, partial derivatives are performed on the odd polynomial expression and the even polynomial expression respectively, and the seventh partial derivative relationship between the odd polynomial and the layer thickness and the eighth partial derivative relationship between the even polynomial and the layer thickness are obtained. The seventh partial derivative relationship is expressed as:
[0118]
[0119] The eighth partial derivative relationship is expressed as:
[0120]
[0121] Substituting the seventh partial derivative relationship formula and the eighth partial derivative relationship formula into the fourth partial derivative calculation formula, the second partial derivative relationship formula of the kernel function for the thickness of each layer of soil can be obtained.
[0122] In a preferred embodiment of the present application, odd polynomials and even polynomials are introduced to represent the Hankel integral kernel function of the first layer, and the odd polynomials and even polynomials are represented by exhaustive permutations and combinations to simplify the partial derivatives of the Hankel integral kernel function of the first layer with respect to soil resistivity and layer thickness. A method for obtaining the partial derivatives of the Hankel integral kernel function of the first layer with respect to soil resistivity and layer thickness is provided, which is easy to program, systematize and implement, and effectively reduces the complexity of the algorithm.
[0123] S3. According to the first partial derivative relationship and the second partial derivative relationship, the first partial derivative calculation formula and the second partial derivative calculation formula are respectively subjected to complex mirror fitting by the complex mirror method, and the third partial derivative relationship and the fourth partial derivative relationship are obtained correspondingly. The first partial derivative calculation formula and the third partial derivative relationship are both used for the partial derivative calculation of the inversion objective function for the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for the partial derivative calculation of the inversion objective function for the layer thickness. In a preferred embodiment of the present application, according to the fifth partial derivative calculation formula and the sixth partial derivative calculation formula, it can be known that the first kernel function correlation relationship related to the kernel function in the fifth partial derivative calculation formula and the second kernel function correlation relationship related to the kernel function in the sixth partial derivative calculation formula are subjected to complex mirror fitting, and the third partial derivative relationship of the inversion objective function for the soil resistivity partial derivative and the fourth partial derivative relationship of the inversion objective function for the layer thickness partial derivative can be further obtained. The first kernel function correlation relationship related to the kernel function in the fifth partial derivative calculation formula is subjected to complex mirror fitting to obtain the first complex mirror fitting result:
[0124]
[0125]
[0126] Perform complex mirror fitting on the second kernel function correlation equation related to the kernel function in the sixth partial derivative calculation formula to obtain the second complex mirror fitting result:
[0127]
[0128] in, is the number of complex images, Indicates the number of the complex image, g is an integer between 1 and n, Express about The number of duplicate images, express The complex image amplitude of express The complex mirror mode of express ( >1) the number of duplicate mirrors, express ( >1), express ( >1) complex mirror mode, express The number of duplicate images, express The complex image amplitude of express The complex mirror mode of .
[0129] Substituting the first complex mirror fitting result into the fifth partial derivative calculation formula, we get:
[0130]
[0131] Substituting the first partial derivative calculation formula further, we can obtain the third partial derivative relationship of the inversion objective function to soil resistivity:
[0132]
[0133] Substituting the second complex mirror fitting result into the sixth partial derivative calculation formula, we get:
[0134]
[0135] Substituting the second partial derivative calculation formula further, we can obtain the fourth partial derivative relationship of the inversion objective function to the layer thickness:
[0136]
[0137] This application fully exports The partial derivative of the layered Hankel integral kernel function on the soil resistivity value and the The partial derivative of the layered Hankel integral kernel function with respect to the layer thickness is used to obtain the partial derivative of the inversion objective function with respect to soil resistivity and the partial derivative of the inversion objective function with respect to layer thickness. The calculation speed and accuracy are improved, ensuring the high efficiency of the inversion of horizontal layered geodetic parameters in the target geodetic area and the accuracy of the inversion results.
[0138] S4. Based on the third partial derivative relationship and the fourth partial derivative relationship, a trust region algorithm is used to perform horizontal layered geodetic parameter inversion to obtain the soil resistivity value and layer thickness value of each layer in the target geodetic region. In a preferred embodiment of the present application, a trust region algorithm is used to perform horizontal layered geodetic parameter inversion. The basic idea of the trust region algorithm is to use a given trust region radius as the upper limit of the displacement length of the soil resistivity and layer thickness, and to determine a spherical region as the trust region region with the current iteration point as the center and the trust region radius. The candidate displacement is determined by solving the optimal point of the quadratic approximation model of the inversion objective function within the trust region region. If the candidate displacement can cause the inversion objective function value to decrease sufficiently, the candidate displacement is accepted as the new displacement, and the trust region radius is maintained or expanded, and a new iteration is continued. Otherwise, it means that the approximation between the quadratic model and the inversion objective function is not ideal, and the trust region radius needs to be reduced. Then, a new candidate displacement is obtained by solving the new optimal point until the iteration termination condition is met.
[0139] Specifically, select the initial parameters 、 、 、 and ,in, and Represent the lower and upper limit proportional coefficients respectively, and denote the contraction coefficient and expansion coefficient respectively, represents the error control amount, , , , in this application, , 、 、 .
[0140] Furthermore, the inversion point vector is constructed based on soil resistivity and soil thickness, and the inversion starting point vector is determined as , , express dimensional real number space.
[0141] Furthermore, the initial gradient vector of the inversion objective function is constructed according to the partial derivative of soil resistivity and the partial derivative of layer thickness corresponding to each layer. The upper limit of the trust region radius and the initial trust region radius are set according to the initial gradient vector of the inversion objective function. The calculation formula is:
[0142]
[0143]
[0144] in, represents the upper limit of the trust region radius, represents the initial gradient vector of the inversion objective function, represents the initial trust region radius, Represents the vector 2-norm.
[0145] According to the third partial derivative relation and the fourth partial derivative relation, the current gradient vector of the inversion objective function is constructed, and the current gradient vector of the inversion objective function is expressed as:
[0146]
[0147] in, represents the number of inversion iterations, Indicates the The current gradient vector of the inversion objective function of the iteration, Indicates the The inversion point vector of the inversion iteration, Represents the gradient.
[0148] Set the convergence judgment condition to When the convergence judgment condition is met, the iteration ends and the output is .
[0149] Specifically, the current gradient vector of the inversion objective function is calculated, and the current gradient vector of the inversion objective function is judged according to the convergence judgment condition. When the current gradient vector of the inversion objective function meets the convergence judgment condition, the output is , and according to The soil resistivity and soil thickness values of each layer in the target earth area are obtained.
[0150] If the current gradient vector of the objective function does not meet the preset convergence judgment condition, the optimal step length is solved within the current trust region radius according to the third and fourth partial derivative relations. The calculation formula for the optimal step length is:
[0151]
[0152] in, Indicates the The optimal step size of the inversion iteration, represents the search range column vector, Indicates the The radius of the dependency domain of the inversion iteration, Represents the layered soil parameters to search for The objective function value.
[0153] According to the optimal step length, the ratio of the actual drop of the inversion objective function to the predicted drop of the inversion objective function is calculated. The calculation formula is:
[0154]
[0155] in, Indicates the The matrix of feasible search directions for the inversion iteration.
[0156] The current trust region radius is iteratively corrected according to the ratio. The iterative correction formula is:
[0157]
[0158] Determine whether the ratio is greater than the preset lower limit proportional coefficient. If so, update the current inversion point. , accept the updated inversion point and update the feasible search direction matrix. The update formula of the feasible search direction matrix is:
[0159]
[0160] Let s = s + 1, update the number of inversion iterations, return to the iterative inversion process and continue the iterative operation until the inversion result is output , and according to The soil resistivity and soil thickness values of each layer in the target earth area are obtained.
[0161] If not, keep the current inversion point, i.e. , set s=s+1, update the number of inversion iterations, return to the iterative inversion process and continue the iterative operation until the inversion result is output , and according to The soil resistivity and soil thickness values of each layer in the target earth area are obtained.
[0162] In a preferred embodiment of the present invention, the apparent resistivity measurement value and the corresponding pole distance of the target earth area are obtained, and the inversion target function is constructed according to the apparent resistivity prediction function and the apparent resistivity measurement value corresponding to the pole distance constructed by using the kernel function and the first kind of zero-order Bessel function; the kernel function is represented by odd polynomials and even polynomials, and a permutation and combination algorithm is used to obtain several ordered sets, and the odd polynomials and even polynomials are represented by several ordered sets and then partial derivatives are calculated, and the first partial derivative relationship of the kernel function to the soil resistivity and the second partial derivative relationship of the kernel function to the layer thickness are obtained, and the ordered sets include the change of soil resistivity of adjacent layers. rate set and layer thickness and set; according to the first partial derivative relationship and the second partial derivative relationship, the complex mirror method is used to perform complex mirror fitting on the first partial derivative calculation formula and the second partial derivative calculation formula respectively, and the third partial derivative relationship and the fourth partial derivative relationship are obtained accordingly. The first partial derivative calculation formula and the third partial derivative relationship are both used to calculate the partial derivative of the inversion target function for soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used to calculate the partial derivative of the inversion target function for layer thickness; according to the third partial derivative relationship and the fourth partial derivative relationship, the trust region algorithm is used to perform horizontal layered earth parameter inversion to obtain the soil resistivity value and layer thickness value of each layer in the target earth area. The horizontal layered geodetic parameter inversion method disclosed in this application uses the Hankel integral kernel function to describe the axisymmetric field problem of radial diffusion of current in the stratum in the four-probe method, converting the complex spatial domain problem into algebraic operations in the wavenumber domain, greatly simplifying the response calculation of multi-layered soil media. The smooth mathematical properties of the first-kind zero-order Bessel function avoid oscillation or divergence problems in numerical calculations, maintain the stability of the apparent resistivity prediction function, and improve the accuracy and reliability of the inversion results. The introduction of odd and even polynomials, combined with the permutation and combination algorithm, effectively improves the calculation speed and accuracy of the partial derivative of the inversion objective function with respect to soil resistivity and the partial derivative of the inversion objective function with respect to layer thickness, ensuring the high efficiency of horizontal layered geodetic parameter inversion in the target geodetic region and the accuracy of the inversion results.
[0163] Accordingly, if Figure 2 FIG2 is a schematic diagram of the structure of a horizontally layered geodetic parameter inversion system. Based on a horizontally layered geodetic parameter inversion method, an embodiment of the present invention further provides a horizontally layered geodetic parameter inversion system to implement the horizontally layered geodetic parameter inversion method disclosed in an embodiment of the present invention, comprising: an inversion objective function construction module 1, a kernel function partial derivative deduction module 2, an inversion objective function partial derivative deduction module 3, and a horizontally layered geodetic parameter inversion module 4.
[0164] The inversion objective function construction module 1 is used to obtain the apparent resistivity measurement value and the corresponding pole distance of the target earth area, and construct the inversion objective function according to the apparent resistivity prediction function corresponding to the pole distance constructed using the kernel function and the first-kind zero-order Bessel function and the apparent resistivity measurement value;
[0165] The kernel function partial derivative deduction module 2 is used to express the kernel function using odd polynomials and even polynomials, and use a permutation and combination algorithm to obtain a plurality of ordered sets. The odd polynomials and the even polynomials are expressed using the plurality of ordered sets and then partial derivatives are calculated to obtain a first partial derivative relationship between the kernel function and soil resistivity and a second partial derivative relationship between the kernel function and layer thickness. The ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums.
[0166] The inversion objective function partial derivative deduction module 3 is used to perform complex mirror fitting on the first partial derivative calculation formula and the second partial derivative calculation formula respectively according to the first partial derivative relationship and the second partial derivative relationship, and obtain a third partial derivative relationship and a fourth partial derivative relationship respectively. The first partial derivative calculation formula and the third partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the layer thickness;
[0167] The horizontal layered earth parameter inversion module 4 is used to perform horizontal layered earth parameter inversion based on the third partial derivative relationship and the fourth partial derivative relationship using a trust region algorithm to obtain the soil resistivity value and layer thickness value of each layer in the target earth area.
[0168] For the specific definition of a horizontally layered earth parameter inversion system, please refer to the above-mentioned definition of a horizontally layered earth parameter inversion method, which will not be repeated here. Those skilled in the art will appreciate that the various modules and steps described in conjunction with the embodiments disclosed in the present invention can be implemented in hardware, software, or a combination of both. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed the scope of the present invention.
[0169] like Figure 3 The internal structure diagram of the computer device shown in the figure, an embodiment of the present invention provides a computer device, including a processor, a memory, and a computer program stored in the memory and configured to be executed by the processor. When the processor executes the computer program, the steps in the embodiment of the horizontal layered earth parameter inversion method described above are implemented, for example Figure 1 Steps S1 to S4 described in .
[0170] Those skilled in the art will understand that the schematic Figure 3 These are merely examples of computer devices and do not constitute limitations on the computer device. The computer device may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the computer device may also include input and output devices, network access devices, buses, etc.
[0171] The processor may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc. The processor is the control center of the computer device, connecting various parts of the entire computer device using various interfaces and lines.
[0172] The memory can be used to store the computer programs and / or modules. The processor implements the various functions of the computer device by running or executing the computer programs and / or modules stored in the memory and accessing the data stored in the memory. The memory may mainly include a program storage area and a data storage area. The program storage area may store an operating system and at least one application required for a function (such as a sound playback function, an image playback function, etc.); the data storage area may store data generated based on the use of the mobile phone (such as audio data, a phone book, etc.). In addition, the memory may include high-speed random access memory and non-volatile memory, such as a hard disk, internal memory, a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, at least one disk storage device, a flash memory device, or other volatile solid-state storage device.
[0173] If the module integrated into the computer device is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the present invention can also implement all or part of the process steps in the above-mentioned method embodiments by instructing the relevant hardware through a computer program. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signal, telecommunication signal, and software distribution medium.
[0174] Those skilled in the art will appreciate that all or part of the processes in the above-described method embodiments can be implemented by instructing the relevant hardware through a computer program. The program can be stored in a computer-readable storage medium, and when executed, the program can include the processes in the above-described method embodiments. The storage medium can be a magnetic disk, an optical disk, a read-only memory (ROM), or a random access memory (RAM).
[0175] Accordingly, an embodiment of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium includes a stored computer program, wherein when the computer program is executed, the device where the computer-readable storage medium is located is controlled to perform the steps in the embodiment of the horizontal layered earth parameter inversion method as described above, for example Figure 1 Steps S1 to S4 described in .
[0176] In summary, the embodiments of the present application provide a method, system, device and medium for inversion of horizontally layered earth parameters, which solve the technical problem of how to improve the accuracy of inversion of horizontally layered earth parameters. The method includes: obtaining the apparent resistivity measurement value and the corresponding pole distance of the target earth area, and constructing an inversion target function according to the apparent resistivity prediction function and the apparent resistivity measurement value corresponding to the pole distance constructed by using the kernel function and the first-class zero-order Bessel function; using odd polynomials and even polynomials to represent the kernel function, using a permutation and combination algorithm to obtain several ordered sets, and using several ordered sets to represent the odd polynomials and even polynomials and then performing partial derivative calculations, and correspondingly obtaining the first partial derivative relationship of the kernel function to the soil resistivity and the kernel function to the partial derivative relationship. The second partial derivative relationship of layer thickness, the ordered set includes the set of adjacent layer soil resistivity change rates and the set of layer thickness sums; according to the first partial derivative relationship and the second partial derivative relationship, the complex mirror method is used to perform complex mirror fitting on the first partial derivative calculation formula and the second partial derivative calculation formula respectively, and the third partial derivative relationship and the fourth partial derivative relationship are obtained accordingly. The first partial derivative calculation formula and the third partial derivative relationship are both used to calculate the partial derivative of the inversion target function for soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used to calculate the partial derivative of the inversion target function for layer thickness; according to the third partial derivative relationship and the fourth partial derivative relationship, the trust region algorithm is used to perform horizontal layered earth parameter inversion to obtain the soil resistivity value and layer thickness value of each layer in the target earth area. The horizontal layered geodetic parameter inversion method disclosed in this application uses the Hankel integral kernel function to describe the axisymmetric field problem of radial diffusion of current in the stratum in the four-probe method, converting the complex spatial domain problem into algebraic operations in the wavenumber domain, greatly simplifying the response calculation of multi-layered soil media. The smooth mathematical properties of the first-kind zero-order Bessel function avoid oscillation or divergence problems in numerical calculations, maintain the stability of the apparent resistivity prediction function, and improve the accuracy and reliability of the inversion results. The introduction of odd and even polynomials, combined with the permutation and combination algorithm, effectively improves the calculation speed and accuracy of the partial derivative of the inversion objective function with respect to soil resistivity and the partial derivative of the inversion objective function with respect to layer thickness, ensuring the high efficiency of horizontal layered geodetic parameter inversion in the target geodetic region and the accuracy of the inversion results.
[0177] Each embodiment in this specification is described in a progressive manner, and the same or similar parts of each embodiment can be referred to each other, and each embodiment focuses on the differences from other embodiments. In particular, for the system embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiment. It should be noted that the various technical features of the above embodiments can be combined arbitrarily. In order to make the description concise, not all possible combinations of the various technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.
[0178] The above-described embodiments merely represent several preferred implementations of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent application. It should be noted that a person skilled in the art could make several improvements and substitutions without departing from the technical principles of the present application, and such improvements and substitutions should also be considered within the scope of protection of the present application. Therefore, the scope of protection of the present patent application shall be based on the scope of protection of the claims.
Claims
1. A horizontal layered earth parameter inversion method, characterized in that: The method comprises: Obtaining apparent resistivity measurement values and corresponding pole distances of a target earth region, and constructing an inversion objective function based on an apparent resistivity prediction function corresponding to the pole distance constructed using a kernel function and a first-kind zero-order Bessel function and the apparent resistivity measurement values; The kernel function is represented by odd polynomials and even polynomials, and a permutation and combination algorithm is used to obtain a plurality of ordered sets. The odd polynomials and the even polynomials are represented by the plurality of ordered sets and then partial derivatives are calculated to obtain a first partial derivative relationship between the kernel function and soil resistivity and a second partial derivative relationship between the kernel function and layer thickness, wherein the ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums; According to the first partial derivative relationship and the second partial derivative relationship, the first partial derivative calculation formula and the second partial derivative calculation formula are respectively subjected to complex mirror fitting using a complex mirror method to obtain a third partial derivative relationship and a fourth partial derivative relationship, wherein the first partial derivative calculation formula and the third partial derivative relationship are both used for calculating the partial derivative of the inversion objective function with respect to the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for calculating the partial derivative of the inversion objective function with respect to the layer thickness; According to the third partial derivative relationship and the fourth partial derivative relationship, a trust region algorithm is used to perform horizontal layered earth parameter inversion to obtain soil resistivity values and layer thickness values of each layer in the target earth region.
2. The horizontal layered earth parameter inversion method according to claim 1, characterized in that: The permutation and combination algorithm is used to obtain a plurality of ordered sets, and the odd polynomial and the even polynomial are represented by the plurality of ordered sets and then partial derivative calculations are performed to obtain a first partial derivative relationship of the kernel function to soil resistivity and a second partial derivative relationship of the kernel function to layer thickness, including: Constructing a first set according to the soil layer number, setting the number of iterations, and using a permutation and combination algorithm to arbitrarily select the number of elements from the first set and arrange them in ascending order to obtain a second set; Using the elements of the second set as numbers, and obtaining the soil resistivity change rate of adjacent layers corresponding to each number according to a preset first calculation rule, permuting and combining the numbers of the soil resistivity change rates of adjacent layers in ascending order to obtain the soil resistivity change rate set of adjacent layers, obtaining the layer thickness sum corresponding to each number according to a preset second calculation rule, permuting and combining the numbers of the layer thickness sums in ascending order to obtain the layer thickness sum set, obtaining the segmented value corresponding to each number according to a preset segmented value assignment rule, and permuting and combining the numbers of the segmented value values in ascending order to obtain a third set; Based on the set of soil resistivity change rates of adjacent layers and the set of layer thickness sums, the odd polynomial is represented according to the properties of the odd polynomial and then iteratively calculated to obtain an odd polynomial expression; and the even polynomial is represented according to the properties of the even polynomial and then iteratively calculated to obtain an even polynomial expression; Obtaining a third partial derivative calculation formula of the kernel function with respect to the soil resistivity, and updating the third partial derivative calculation formula according to the odd polynomial expression and the even polynomial expression to obtain a first partial derivative relationship formula of the kernel function with respect to the soil resistivity; The fourth partial derivative calculation formula of the kernel function for the layer thickness is updated according to the odd polynomial expression, the even polynomial expression and the third set to obtain the second partial derivative relationship of the kernel function for the layer thickness.
3. The horizontal layered earth parameter inversion method according to claim 2, characterized in that: The odd polynomial expression is: The even polynomial expression is: in, represents the integration variable, represents the number of iterations, Indicates the total number of soil stratification layers, represents the elements of the second set, represents the first element of the second set, represents the second element of the second set, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the numbered layer thickness and, represents the set of soil resistivity change rates of adjacent layers, represents the layer thickness and set, Indicates from Choose from natural numbers The number of combinations of natural numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the change rate of soil resistivity of adjacent layers with numbers, Represents the elements of the second set is the numbered layer thickness and, Represents the elements of the second set is the thickness and of the numbered layers.
4. The horizontal layered earth parameter inversion method according to claim 1, wherein: The apparent resistivity prediction function is expressed as: in, Indicates the The measured pole distance, Indicates pole distance The corresponding predicted apparent resistivity value is, Indicates the Layered soil resistivity, Indicates the Hierarchical Hankel integral kernel function, express The corresponding zero-order Bessel function of the first kind is, express The corresponding zero-order Bessel function of the first kind is, represents the integration variable.
5. The horizontal layered earth parameter inversion method according to claim 1, characterized in that: The inversion objective function is constructed based on the apparent resistivity prediction function corresponding to the pole distance constructed using a kernel function and a first-kind zero-order Bessel function and the apparent resistivity measurement value, including: Obtaining a measurement error according to the apparent resistivity prediction function and the apparent resistivity measurement value; An inversion objective function is constructed with the goal of minimizing the ratio of the measurement error to the apparent resistivity measurement value.
6. The horizontal layered earth parameter inversion method according to claim 1, wherein: According to the first partial derivative relationship and the second partial derivative relationship, the first partial derivative calculation formula and the second partial derivative calculation formula are respectively subjected to complex mirror fitting by the complex mirror method, and a third partial derivative relationship and a fourth partial derivative relationship are correspondingly obtained, including: Obtaining a first kernel function correlation formula in a fifth partial derivative calculation formula of the apparent resistivity prediction function with respect to the soil resistivity, and a second kernel function correlation formula in a sixth partial derivative calculation formula of the apparent resistivity prediction function with respect to the layer thickness; According to the first partial derivative relationship and the second partial derivative relationship, a complex mirror method is used to perform complex mirror fitting on the first kernel function correlation relationship and the second kernel function correlation relationship, respectively, to obtain a first complex mirror fitting result and a second complex mirror fitting result; Obtaining a first partial derivative calculation formula of the inversion objective function with respect to the soil resistivity, and sequentially updating the fifth partial derivative calculation formula and the first partial derivative calculation formula according to the first complex mirror fitting result to obtain a third partial derivative relationship formula of the inversion objective function with respect to the soil resistivity; Obtain a second partial derivative calculation formula of the inversion objective function for the layer thickness, and update the sixth partial derivative calculation formula and the second partial derivative calculation formula in sequence according to the second complex mirror fitting result to obtain a fourth partial derivative relationship of the inversion objective function for the layer thickness.
7. The horizontal layered earth parameter inversion method according to claim 1, characterized in that: The method of performing horizontal layered earth parameter inversion based on the third partial derivative relationship and the fourth partial derivative relationship using a trust region algorithm to obtain soil resistivity values and layer thickness values of each layer in the target earth region includes: Determine the inversion starting point, and obtain the soil resistivity partial derivative of the inversion objective function for each layer according to the third partial derivative relationship, and obtain the layer thickness partial derivative of the inversion objective function for each layer according to the fourth partial derivative relationship; Obtaining an initial gradient vector of an inversion objective function according to the partial derivative of the soil resistivity and the partial derivative of the layer thickness corresponding to each layer, and setting an upper limit of a trust region radius and an initial trust region radius according to the initial gradient vector of the objective function; During the iterative inversion process, a current gradient vector of the inversion objective function is calculated, and it is determined whether the current gradient vector of the inversion objective function satisfies a preset convergence judgment condition. If so, an inversion result is output, wherein the inversion result includes soil resistivity values and soil thickness values of each layer of the target earth region; If not, solving the optimal step length within the current trust region radius limit according to the third partial derivative relationship and the fourth partial derivative relationship, and calculating the ratio of the actual descent amount of the inversion objective function to the predicted descent amount of the inversion objective function according to the optimal step length, and iteratively correcting the current trust region radius according to the ratio; Determining whether the ratio is greater than a preset lower limit proportional coefficient, if so, updating the current inversion point and updating the feasible search direction matrix based on the updated current inversion point, if not, maintaining the current inversion point; The number of iterations is updated to return to the iterative inversion process to continue iterative operations until the inversion result is output.
8. A horizontal layered earth parameter inversion system, used to implement the horizontal layered earth parameter inversion method according to any one of claims 1 to 7, characterized in that: The system comprises: an inversion objective function construction module, a kernel function partial derivative deduction module, an inversion objective function partial derivative deduction module and a horizontal layered earth parameter inversion module; The inversion objective function construction module is used to obtain the apparent resistivity measurement value and the corresponding pole distance of the target earth area, and construct the inversion objective function according to the apparent resistivity prediction function corresponding to the pole distance constructed using the kernel function and the first-kind zero-order Bessel function and the apparent resistivity measurement value; The kernel function partial derivative deduction module is used to represent the kernel function using odd polynomials and even polynomials, and to obtain a plurality of ordered sets using a permutation and combination algorithm. The odd polynomials and the even polynomials are represented using the plurality of ordered sets and then partial derivative calculations are performed to obtain a first partial derivative relationship between the kernel function and soil resistivity and a second partial derivative relationship between the kernel function and layer thickness. The ordered sets include a set of soil resistivity change rates of adjacent layers and a set of layer thickness sums. The inversion objective function partial derivative deduction module is used to perform complex mirror fitting on the first partial derivative calculation formula and the second partial derivative calculation formula respectively according to the first partial derivative relationship and the second partial derivative relationship using the complex mirror method, and obtain a third partial derivative relationship and a fourth partial derivative relationship respectively. The first partial derivative calculation formula and the third partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the soil resistivity, and the second partial derivative calculation formula and the fourth partial derivative relationship are both used for calculating the partial derivative of the inversion objective function for the layer thickness; The horizontal layered earth parameter inversion module is used to perform horizontal layered earth parameter inversion based on the third partial derivative relationship and the fourth partial derivative relationship using a trust region algorithm to obtain the soil resistivity value and layer thickness value of each layer in the target earth area.
9. A computer device, characterized in that: The computer device includes a memory, a processor and a transceiver, which are connected via a bus; the memory is used to store a set of computer program instructions and data, and transmit the stored data to the processor, and the processor executes the computer program instructions stored in the memory to perform the horizontal layered earth parameter inversion method as described in any one of claims 1 to 7.
10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program, and when the computer program is executed, the horizontal layered earth parameter inversion method according to any one of claims 1 to 7 is implemented.
Citation Information
Patent Citations
Magnetotelluric regularization inversion method based on different constraint conditions
CN104360404A
Earth resistivity model modeling method and device, computer equipment and storage medium
CN114047554A