A fast calculation method for remote detection of complex formation structure while drilling electromagnetic wave logging response

Abstract

The present invention discloses a method for rapidly calculating the response of remotely detected while-drilling electromagnetic wave logging in complex formation structures, relating to the field of electrical logging technology for oil exploration and development. The method comprises the following steps: establishing a two-dimensional heterogeneous isotropic geological model of the complex formation structure; using a regional rapid and precise description method to divide the formation into a collection of triangular units and assign resistivity values; determining the calculation region through adaptive thresholding, segmenting the x and z directions, and sequentially assigning resistivity values ​​to the calculation region grids according to the segmentation points to generate a calculation model; obtaining a large sparse matrix linear equation system; solving the equation system using a hybrid solver to obtain the distribution of the spectral domain field; calculating the integral of the spectral domain field with wave number according to the inverse Fourier transform formula to obtain the spatial domain distribution of the electromagnetic field; and moving each logging point to obtain the while-drilling electromagnetic wave logging response curve for the complex formation environment. The present invention can effectively and rapidly calculate the electromagnetic field response of complex two-dimensional formations, improving the algorithm's adaptability and calculation speed.

Description

Technical Field

[0001] The present invention relates to the technical field of electrical logging for oil exploration and development, and belongs to the category of numerical forward calculation methods, in particular to a fast calculation method for remote detection while-drilling electromagnetic wave logging responses for complex stratum structures. Background Art

[0002] As exploration depths continue to increase and the demand for azimuth signals grows, one-dimensional forward modeling is increasingly unable to meet the demands of long-distance exploration. Three-dimensional forward modeling lacks analytical solutions, and numerical calculations are slow, making it impossible to meet production requirements. Since existing instruments typically have a detection range of around 30 meters, the formation can be simplified into a two-dimensional layer with electrical properties constant in one direction and varying only in the other two. Therefore, 2.5-dimensional numerical forward modeling can better meet the computational speed and accuracy requirements for simulating the distribution of three-dimensional electromagnetic fields.

[0003] In order to improve the solution speed of sparse stiffness matrix and adapt to more complex two-dimensional formations, it is urgent to study a fast calculation method for remote detection of while-drilling electromagnetic wave logging for complex formation structures. Summary of the Invention

[0004] In order to solve the above technical problems, the present invention discloses a method for rapid calculation of the response of electromagnetic wave logging while drilling for remote detection of complex formation structures. The method adopts a regional rapid and accurate description method to construct complex two-dimensional formations, and accelerates the iterative solution through methods such as adaptive segmentation boundary truncation and sparse matrix hybrid solution algorithm, so as to quickly solve the electromagnetic field distribution through the frequency domain finite difference algorithm, accelerate the numerical forward modeling speed of two-dimensional complex formations, and provide strong support for remote detection of electromagnetic wave logging while drilling.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A method for rapidly calculating the response of electromagnetic wave logging while drilling for remote detection of complex formation structures comprises the following steps:

[0007] s1. Establish a two-dimensional heterogeneous isotropic geological model of complex stratigraphic structures, set the wellbore trajectory, determine the position coordinates of each logging point in the trajectory, and set instrument parameters;

[0008] s2. Use the regional fast and accurate description method to divide the stratum into a set of triangular units and assign a resistivity value to each triangular unit;

[0009] s3. Calculate the field distribution characteristics based on the plane wave attenuation law, determine the calculation area through threshold adaptation, partition the optimized calculation area in the x and z directions, assign resistivity values ​​to the calculation area grids according to the partition points, and finally generate the calculation model;

[0010] s4. Convert the spatial Helmholtz equations to the wavenumber domain using Fourier transform and install them using the Yee cellular finite difference scheme to obtain the large sparse matrix linear equations to be solved;

[0011] s5. Solve large sparse matrix linear equations using a hybrid solver, quickly and accurately calculating the distribution of the spectral field;

[0012] s6. Calculate the integral of the spectral field with wave number according to the inverse Fourier transform formula to obtain the spatial distribution of the electromagnetic field;

[0013] s7. Move each logging point and loop through steps s3 to s6 to obtain a while-drilling electromagnetic wave logging response curve in a complex formation environment.

[0014] Optionally, in step s2, the method for quickly and accurately describing the region is:

[0015] s21. Construct a two-dimensional plane range, set boundaries, determine the number of regions in the two-dimensional plane, and set the number of points, lines, and planes. Set the position coordinates of the points, the points corresponding to each line, and the lines contained in each plane.

[0016] s22. Divide the surface area into triangular areas and divide the two-dimensional plane by determining the position of the points and the triangular areas, and set the resistivity of different triangular areas;

[0017] s23. When assigning grid resistivity, each grid is divided into smaller cellular structures. The position of each cell in the triangular area is searched point by point in the cellular order to obtain the cellular resistivity.

[0018] s24. Perform equivalent processing on the cell resistivity values ​​obtained in step s23 to obtain the resistivity of the corresponding grid, and form a grid resistivity matrix after processing all grids.

[0019] Optionally, in step s3, the calculation area is determined by threshold adaptively:

[0020] s31. First, determine the maximum calculation range by calculating the instrument source distance, frequency, and resistivity of the launch point. Use a coarse gridding method that covers the entire area to construct the initial calculation area.

[0021] s32. Calculate the electromagnetic wave attenuation coefficient in each grid based on the propagation characteristics of electromagnetic waves. Multiply the electromagnetic field attenuation coefficients of consecutive grid cells in the horizontal direction to calculate the electromagnetic wave attenuation value at each grid node.

[0022] s33. Using the electromagnetic wave attenuation value at the transverse grid node as the base value, multiply the plane wave attenuation value along the longitudinal grid and compare it with the threshold. If it does not reach the threshold, continue the multiplication; if it exceeds the threshold, terminate the multiplication.

[0023] s34. Monitor the relationship between the accumulated value and the field strength threshold in real time, and record the corresponding extreme value of the longitudinal propagation distance when the condition is met. After traversing all longitudinal grids, compare the extreme values ​​corresponding to each benchmark parameter, and select the global maximum value as the basis for determining the effective propagation range of the electromagnetic wave and the boundary distance of the calculation area;

[0024] s35. Based on the newly determined calculation area boundary, regenerate the grid to generate a calculation grid that adapts to the formation model.

[0025] Optionally, in step s5, a large sparse matrix linear equation system is solved by a hybrid solver, specifically:

[0026] s51. Set the wavenumber sampling array according to the spectral field convergence characteristics, and set the wavenumber threshold k according to the calculation speed of the direct solver and the iterative solver by comparing their computational efficiency. y0 ;

[0027] s52. Determine the size of the current wave value and the wave number cutoff point. If it is less than the wave number threshold k y0 , then solve it by direct solver;

[0028] s53. If the current wave value is greater than the wave number threshold k y0 , then it is solved by iterative solver;

[0029] s54. According to the given wavenumber sampling array, for each wavenumber point in the array, loop through steps s52 and s53 to perform judgment and solution, and obtain all spectral domain field solutions in the entire spectral domain.

[0030] The beneficial effect of the present invention is that the present invention provides a method for quickly calculating the response of electromagnetic wave logging while drilling for remote detection of complex formation structures. The two-dimensional complex formation description method, the adaptive truncation of the partitioning boundary and the sparse matrix hybrid solution method are applied to the traditional 2.5-dimensional finite difference calculation method. The electromagnetic field response of complex two-dimensional formations can be effectively and quickly calculated, which greatly improves the adaptability and calculation speed of the algorithm and provides strong support for remote detection of electromagnetic wave logging. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 This is a flow chart of a method for quickly calculating the response of electromagnetic wave logging while drilling for remote detection of complex formation structures according to the present invention;

[0032] Figure 2 A diagram of a complex two-dimensional geological body model according to an embodiment of the present invention;

[0033] Figure 3 A method for dividing a region into triangular elements is shown as an embodiment of the present invention;

[0034] Figure 4 A schematic diagram of wellbore trajectory selection and calculation area cutting according to an embodiment of the present invention;

[0035] Figure 5 A comparison diagram of mesh subdivision before and after adaptive regional truncation according to an embodiment of the present invention;

[0036] Figure 6 A comparison chart of the solution time of a hybrid solver and other solvers shown in one embodiment of the present invention;

[0037] Figure 7 A diagram showing a rapid calculation response of an electric field original signal according to an embodiment of the present invention;

[0038] Figure 8 A diagram illustrating a rapid calculation response of a magnetic field original signal according to an embodiment of the present invention. DETAILED DESCRIPTION

[0039] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. 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. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention for which protection is sought, but merely represents selected embodiments 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.

[0040] A fast calculation method for remote detection of complex formation structure while drilling electromagnetic wave logging response, such as Figure 1 As shown, the following steps are included:

[0041] s1. Construct a two-dimensional heterogeneous isotropic electrical logging calculation model for complex formation structures (such as faults, blocks, anomalies, etc.) and drilling environments, set the formation resistivity value Rt and the wellbore trajectory l, and determine the position coordinates P of each logging point in the trajectory i (x i ,y i ); Set instrument parameters such as source distance TR, frequency f, etc.

[0042] s2. Use a regional fast and accurate description method to divide the stratum into a set of triangular units and assign a resistivity value to each triangular unit. The regional fast and accurate description method is as follows:

[0043] s21. Construct a two-dimensional plane range and set the boundary to (x min ,x max ) and (z min ,z max ), determine the number of regions N in the two-dimensional plane k , set the coordinates of the point P according to the number of regions i (x i ,z i ), line segment l i Points and surfaces S connected to it i The included lines represent the various areas;

[0044] s22. Divide the surface area into triangular areas and divide the two-dimensional plane by determining the position of the points and the triangular areas, and set the resistivity of different triangular areas;

[0045] s23. When assigning grid resistivity, each grid is divided into smaller cellular structures. The position of each cell in the triangular area is searched point by point in the cellular order to obtain the cellular resistivity.

[0046] s24. Perform equivalent processing on the cell resistivity values ​​obtained in step s23 to obtain the resistivity of the corresponding grid, and form a grid resistivity matrix after processing all grids.

[0047] s3. Calculate the distribution characteristics of electromagnetic field attenuation by opening a window based on the instrument inclination angle φ and different background formation resistivities Rt. Calculate the field distribution characteristics based on the plane wave attenuation law, determine the calculation area through threshold adaptation, and then partition the optimized calculation area in the x and z directions. Assign resistivity values ​​to the calculation area grids according to the partition points, and finally generate a calculation model. Specifically, the calculation area determined by threshold adaptation is:

[0048] s31. First, determine the maximum calculation range by calculating the instrument source distance, frequency, and resistivity of the launch point. Use a coarse gridding method that covers the entire area to construct the initial calculation area.

[0049] s32. Calculate the electromagnetic wave attenuation coefficient in each grid based on the propagation characteristics of the electromagnetic wave and , perform electromagnetic field attenuation coefficient multiplication on continuous grid units in the horizontal direction, and calculate the electromagnetic wave attenuation value at each grid node ,in,{ } is the benchmark dataset of longitudinal attenuation, is the attenuation coefficient at the kth calculation grid in the horizontal direction;

[0050] s33. Taking the electromagnetic wave attenuation value at the transverse grid node as the base value, the attenuation value of the plane wave along the longitudinal grid is multiplied. ,in, is the attenuation coefficient at the (i, j)th grid cell, is the addition coefficient of the mth vertical grid at the ith horizontal grid distance. Compare with the threshold, if it is less than the threshold, continue the accumulation, if it exceeds the threshold, stop the accumulation;

[0051] s34. Monitor the relationship between the accumulated value and the field strength threshold in real time, and record the corresponding extreme value of the longitudinal propagation distance when the condition is met. After traversing all longitudinal grids, compare the extreme values ​​corresponding to each benchmark parameter, and select the global maximum value as the basis for determining the effective propagation range of the electromagnetic wave and the boundary distance of the calculation area;

[0052] s35. Based on the newly determined calculation area boundary, regenerate the grid to generate a calculation grid that adapts to the formation model.

[0053] s4. Convert the spatial Helmholtz equations to the wavenumber domain using Fourier transform and install them using the Yee cellular finite difference scheme to obtain the large sparse matrix linear equations to be solved;

[0054] s5. Solve large sparse matrix linear equations using a hybrid solution approach, quickly and accurately calculating the distribution of the spectral field; specifically:

[0055] s51. Set the wavenumber sampling array according to the spectral field convergence characteristics, and set the wavenumber threshold k according to the calculation speed of the direct solver and the iterative solver by comparing their computational efficiency. y0 ;

[0056] s52. Determine the size of the current wave value and the wave number cutoff point. If it is less than the wave number threshold k y0 , then solve it by direct solver;

[0057] s53. If the current wave value is greater than the wave number threshold k y0 , then it is solved by iterative solver;

[0058] s54. According to the given wavenumber sampling array, for each wavenumber point in the array, loop through steps s52 and s53 to perform judgment and solution, and obtain all spectral domain field solutions in the entire spectral domain.

[0059] s6. Accumulate the spectral domain electromagnetic fields obtained by solving for different wave numbers, and obtain the spatial domain solution by inverse Fourier transform of the spectral domain solution of the electromagnetic field;

[0060] s7. Move each logging point P i (x i ,y i ), and looping steps s3 to s6, the electromagnetic field response curve of electromagnetic wave logging while drilling in complex formation environment.

[0061] Application Examples

[0062] A fast calculation method for remote detection of complex formation structure while drilling electromagnetic wave logging response, combined with Figures 2 to 8 , the specific steps are as follows:

[0063] s1.Build as Figure 2 The isotropic model of a layered, isolated anomaly is shown. The resistivity of the upper layer is set to 100 Ω.m, the resistivity of the lower layer is set to 10 Ω.m, and the resistivity of the isolated anomaly is set to 1000 Ω.m. The wellbore trajectory is a vertical wellbore, with horizontal distances from the anomaly at L (0.5 m, 1 m, 1.5 m). The instrument source spacing is set to 5 m, and the frequency is set to 100 kHz.

[0064] s2. For the embodiment model, the stratum is divided into a set of 12 triangular units using a regional fast and accurate description method, such as Figure 3 As shown. A resistivity value is assigned to each triangular element; the method for fast and accurate regional description is:

[0065] s21. Construct a two-dimensional plane range and set the boundary to (x min ,x max ) and (z min ,z max ), determine that there are 12 regions in the two-dimensional plane, and set the coordinates of the point P according to the number of regions i (x i ,z i ),i=1,2,…,10, the line segment contains P1-P2,P1-P3,P1-P4,P2-P4,P3-P4,P3-P5,P3-P8,P3-P 10 ,P4-P6,P4-P7,P4-P 10 ,P5-P6,P5-P8,P5-P9,P6-P7,P6-P9,P7-P9,P7-P 10 ,P8-P9,P8-P 10 ,P9-P 10 , including P1P2P4P3, P3P4P6P5, P7P9P8P 10 , to indicate each area;

[0066] s22. Divide the surface area into triangular areas and divide the two-dimensional plane by determining the position of the points and the triangular areas, and set the resistivity of different triangular areas;

[0067] s23. When assigning grid resistivity, each grid is divided into smaller cellular structures. The position of each cell in the triangular area is searched point by point in the cellular order to obtain the cellular resistivity.

[0068] s24. Perform equivalent processing on the cell resistivity values ​​obtained in step s23 to obtain the resistivity of the corresponding grid, and form a grid resistivity matrix after processing all grids.

[0069] s3. Calculate the distribution characteristics of electromagnetic field attenuation by opening windows based on different background formation resistivity Rt. Calculate the field distribution characteristics based on the plane wave attenuation law and determine the calculation area through threshold adaptation. The calculation area formed after the model adaptive windowing is as follows: Figure 4 As shown, the x and z directions are divided in the optimized calculation area, and the resistivity of the model grid is assigned in sequence according to the division points, and finally the calculation model is generated. The grid division of the calculation area before and after optimization is as follows Figure 5 As shown in the figure, the optimized window boundary is 11.3 m; the calculation area is determined by threshold adaptation as follows:

[0070] s31. First, determine the maximum calculation range by calculating the instrument source distance, frequency, and resistivity of the launch point. Use a coarse gridding method that covers the entire area to construct the initial calculation area.

[0071] s32. Calculate the electromagnetic wave attenuation coefficient in each grid based on the propagation characteristics of the electromagnetic wave and , perform electromagnetic field attenuation coefficient multiplication on continuous grid units in the horizontal direction, and calculate the electromagnetic wave attenuation value at each grid node ,in,{ } is the benchmark dataset of longitudinal attenuation, is the attenuation coefficient at the kth calculation grid in the horizontal direction;

[0072] s33. Taking the electromagnetic wave attenuation value at the transverse grid node as the base value, the attenuation value of the plane wave along the longitudinal grid is multiplied. ,in, is the attenuation coefficient at the (i, j)th grid cell, is the addition coefficient of the mth vertical grid at the ith horizontal grid distance. Compare with the threshold, if it is less than the threshold, continue the accumulation, if it exceeds the threshold, stop the accumulation;

[0073] s34. Monitor the relationship between the accumulated value and the field strength threshold in real time, and record the corresponding extreme value of the longitudinal propagation distance when the condition is met. After traversing all longitudinal grids, compare the extreme values ​​corresponding to each benchmark parameter, and select the global maximum value as the basis for determining the effective propagation range of the electromagnetic wave and the boundary distance of the calculation area;

[0074] s35. Based on the newly determined calculation area boundary, regenerate the grid to generate a calculation grid that adapts to the formation model.

[0075] s4. The spatial Helmholtz equations are converted to the wavenumber domain according to the Fourier transform and installed according to the Yee cellular finite difference scheme to obtain the large sparse matrix linear equations to be solved.

[0076] s5. Solve large sparse matrix linear equations using a hybrid solution approach to quickly and accurately calculate the distribution of the spectral field; specifically:

[0077] s51. Set the wave number sampling array according to the spectral field convergence characteristics, and compare the computational efficiency of the direct solver and the iterative solver, such as Figure 6 As shown, the wave number threshold k is set based on the calculated speed of the two y0 , here k y0 It was obtained at the 11th sampling point, which was about 1.8;

[0078] s52. Determine the size of the current wave value and the wave number cutoff point. If it is less than the wave number threshold k y0 , then solve it by direct solver;

[0079] s53. If the current wave value is greater than the wave number threshold k y0 , then it is solved by iterative solver;

[0080] s54. According to the given wavenumber sampling array, for each wavenumber point in the array, loop through steps s52 and s53 to perform judgment and solution, and obtain all spectral domain field solutions in the entire spectral domain.

[0081] s6. Accumulate the spectral domain electromagnetic fields obtained by solving for different wave numbers, and obtain the spatial domain solution by inverse Fourier transform of the spectral domain solution of the electromagnetic field;

[0082] s7. Move each logging point P i (x i ,y i ), and loop through steps s3 to s6 to obtain Figure 7 and Figure 8The electromagnetic field response curve of electromagnetic wave logging while drilling in complex formation environments shows that this method can effectively simulate the electromagnetic field response of complex geological models such as isolated anomalies. As the distance between the anomaly and the wellbore trajectory changes, the logging response curve changes. The closer the distance between the anomaly and the wellbore, the greater the impact; the farther the distance, the smaller the impact.

[0083] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. A method for rapid calculation of electromagnetic wave logging response while drilling for remote detection of complex formation structures, characterized in that: The steps include: s1. Establish a two-dimensional heterogeneous isotropic geological model of complex stratigraphic structures, set the wellbore trajectory, determine the position coordinates of each logging point in the trajectory, and set instrument parameters; s2. Use the regional fast and accurate description method to divide the stratum into a set of triangular units and assign a resistivity value to each triangular unit; s3. Calculate the field distribution characteristics based on the plane wave attenuation law, determine the calculation area through threshold adaptation, partition the optimized calculation area in the x and z directions, assign resistivity values ​​to the calculation area grids according to the partition points, and finally generate the calculation model; s4. Convert the spatial Helmholtz equations to the wavenumber domain using Fourier transform and install them using the Yee cellular finite difference scheme to obtain the large sparse matrix linear equations to be solved; s5. Solve large sparse matrix linear equations using a hybrid solver to obtain the distribution of the spectral field; s6. Calculate the integral of the spectral field with wave number according to the inverse Fourier transform formula to obtain the spatial distribution of the electromagnetic field; s7. Move each logging point and loop through steps s3 to s6 to obtain a while-drilling electromagnetic wave logging response curve in a complex formation environment.

2. The method for rapid calculation of electromagnetic wave logging response while drilling for remote detection of complex strata according to claim 1, characterized in that: In step s2, the method for quickly and accurately describing the region is: s21. Construct a two-dimensional plane range, set boundaries, determine the number of regions in the two-dimensional plane, and set the number of points, lines, and planes. Set the position coordinates of the points, the points corresponding to each line, and the lines contained in each plane. s22. Divide the surface area into triangular areas and divide the two-dimensional plane by determining the position of the points and the triangular areas, and set the resistivity of different triangular areas; s23. When assigning grid resistivity, divide each grid into a cellular structure, and find its position in the triangular area point by point in the cellular order to obtain the cellular resistivity; s24. Perform equivalent processing on the cell resistivity values ​​obtained in step s23 to obtain the resistivity of the corresponding grid, and form a grid resistivity matrix after processing all grids.

3. The method for rapid calculation of electromagnetic wave logging response while drilling for remote detection of complex strata according to claim 1 is characterized in that: In step s3, the calculation area is determined by threshold adaptively: s31. First, determine the maximum calculation range by calculating the instrument source distance, frequency, and resistivity of the launch point. Use a coarse gridding method that covers the entire area to construct the initial calculation area. s32. Calculate the electromagnetic wave attenuation coefficient in each grid based on the propagation characteristics of electromagnetic waves. Multiply the electromagnetic field attenuation coefficients of consecutive grid cells in the horizontal direction to calculate the electromagnetic wave attenuation value at each grid node. s33. Using the electromagnetic wave attenuation value at the transverse grid node as the base value, multiply the plane wave attenuation value along the longitudinal grid and compare it with the threshold. If it does not reach the threshold, continue the multiplication; if it exceeds the threshold, terminate the multiplication. s34. Monitor the relationship between the accumulated value and the field strength threshold in real time, and record the corresponding extreme longitudinal propagation distance when the condition is met. After traversing all longitudinal grids, compare the extreme values ​​corresponding to each benchmark parameter, and select the global maximum value as the basis for determining the effective electromagnetic wave propagation range and the boundary distance of the calculation area; s35. Based on the newly determined calculation area boundary, regenerate the grid to generate a calculation grid that adapts to the formation model.

4. The method for rapid calculation of electromagnetic wave logging response while drilling for remote detection of complex strata according to claim 1, characterized in that: In step s5, a large sparse matrix linear equation system is solved by a hybrid solver, specifically: s51. Set the wavenumber sampling array according to the spectral field convergence characteristics, and set the wavenumber threshold k according to the calculation speed of the direct solver and the iterative solver by comparing their computational efficiency. y0 ; s52. Determine the size of the current wave value and the wave number cutoff point. If it is less than the wave number threshold k y0 , then solve it by direct solver; s53. If the current wave value is greater than the wave number threshold k y0 , then it is solved by iterative solver; s54. According to the given wavenumber sampling array, for each wavenumber point in the array, loop through steps s52 and s53 to perform judgment and solution, and obtain all spectral domain field solutions in the entire spectral domain.

Images