Electromagnetic inversion method and device based on well electric shock joint modeling

Abstract

The invention relates to geophysical inversion, in particular to an electromagnetic inversion method based on well electric-seismic joint modeling, which comprises the following steps: acquiring actually measured apparent resistivity data; well-to-seismic joint stratum thickness and stratum apparent resistivity modeling based on geosteering; performing one-dimensional band source apparent resistivity inversion based on a nonlinear least square fitting algorithm; and outputting a stratum apparent resistivity inversion result. The method is based on geosteering, the advantages of three geophysical prospecting methods of earthquake, logging and electromagnetism are fully fused and played, the multiplicity of electromagnetic inversion solutions is reduced, the inversion depth and accuracy are improved, and powerful guarantee is provided for subsequent section interpretation, well position deployment, fluid identification and other operations. The invention also provides a device for realizing the electromagnetic inversion method based on well electric shock joint modeling.

Description

Technical Field

[0001] The invention relates to geophysical inversion, and in particular to an electromagnetic inversion method and device based on well-electric-seismic joint modeling. Background Art

[0002] Inversion has always been a difficult problem that needs to be faced in the field of geophysics. Its accuracy is directly related to the smooth progress of subsequent engineering operations such as tunnel hidden danger inspection and oil and gas exploration and development. Electromagnetic exploration and seismic exploration are two important branches of geophysical exploration. Compared with seismic exploration, electromagnetic exploration obtains the apparent resistivity properties corresponding to different frequency points of the detection point. There are fewer frequency points in the vertical direction, and its resolution will drop sharply with the increase of depth. At present, most conventional electromagnetic inversions are free inversions (unconstrained inversions). The inversion depth is determined according to the skin depth calculation formula. It is usually more accurate at a depth of 1 to 2 km. It is mainly used for tunnel detection of roads, railways, and large-scale structural identification such as faults and broken zones. If it is to be applied to oil and gas fields shallower than 5 km, it is difficult to achieve it only by unconstrained inversion due to the limitation of insufficient vertical resolution of single electromagnetic inversion. At the same time, the multi-solution nature of geophysical inversion is also a difficult problem that needs to be faced. Due to the existence of multiple unknowns, hundreds, thousands, or even tens of thousands of different solutions are usually obtained, which seriously affects the subsequent data interpretation.

[0003] In view of this, there is an urgent need to provide a new electromagnetic inversion method that can achieve more efficient and reliable formation apparent resistivity inversion in deeper underground spaces, such as less than 5 km. Summary of the invention

[0004] The purpose of the present invention is to break through the limitations of traditional electromagnetic unconstrained inversion, take seismic spatial structure as the basic framework, and improve the framework in combination with geological knowledge, use logging data to calibrate the apparent resistivity of the formation and set threshold constraints, carry out joint well-electro-seismic modeling inversion, reduce the multi-solution of inversion, and greatly improve the inversion depth and accuracy.

[0005] In order to achieve the above-mentioned object of the invention, the present invention provides an electromagnetic inversion method based on wellbore electric and seismic joint modeling, comprising:

[0006] S1, obtain measured apparent resistivity data;

[0007] S2, well-seismic combined formation thickness and formation apparent resistivity modeling based on geosteering;

[0008] S3, one-dimensional band-source apparent resistivity inversion based on nonlinear least squares fitting algorithm;

[0009] S4, output the formation apparent resistivity inversion results.

[0010] This electromagnetic inversion method can be applied to both traditional electromagnetic exploration methods and wide-area electromagnetic exploration methods. The formation apparent resistivity parameter inversion based on well-electric-seismic joint modeling can significantly improve the effective depth, efficiency and accuracy of electromagnetic inversion. Since the wide-area electromagnetic method is based on the precise solution defined by it, it has the characteristic of a greater effective exploration depth than the traditional electromagnetic method. Taking the wide-area electromagnetic method exploration as an example to illustrate the electromagnetic inversion method based on well-electric-seismic joint modeling, the advantages of the electromagnetic inversion method based on well-electric-seismic joint modeling of the present invention can be better reflected. However, the application scope of the constrained electromagnetic inversion method based on fast forward iteration provided by the present invention is not limited to the wide-area electromagnetic method.

[0011] The measured apparent resistivity data obtained by S1 is the apparent resistivity data calculated by wide-area electromagnetic method, further including:

[0012] S11, performing potential difference preprocessing on the original potential difference data collected by the wide-area electromagnetic method receiver;

[0013] S12, preprocessing the current value of the transmission current at the non-smooth current position of the wide-area electromagnetic method transmitter;

[0014] S13, the measured apparent resistivity value is calculated using the inverse interpolation method.

[0015] The potential difference preprocessing described in S11 includes correcting the abnormal points or areas of the original potential difference by selecting one or any combination of the following methods according to the received preprocessing instructions: median filtering, mean filtering, mean regression, curve translation, and jump point correction. The preprocessing instructions can be the point number, line number, or combination thereof to be processed input by the graphical user interface GUI defined by the computer software program, as well as the interference information of the specific original data acquisition site, and / or the preprocessing method selection input after the abnormal type is identified by manual reading of the map. In other words, for the abnormal points / areas to be preprocessed, the preprocessing method (combination) can be selected based on human experience, or the corresponding interference information-preprocessing method (combination) rules can be preset to realize automatic preprocessing of the potential difference raw data.

[0016] The current value preprocessing described in S12 includes smoothing the connection between different frequency groups of a single measuring point or other non-smooth emission current value abnormal points or areas according to the received preprocessing instructions, so that the frequency group connection and the complete current curve are smooth and continuous; and the complete current curve of the single measuring point is extended to the remaining measuring points of the same field source. Similar to the preprocessing instructions at the original potential difference, the preprocessing instructions here can be the emission current curve data of the point number to be processed input by the graphical user interface GUI defined by the computer software program, and / or the preprocessing method selection input after the current value abnormality type is identified by manual reading of the map.

[0017] After calculating the measured apparent resistivity, S13 also performs apparent resistivity processing, which includes selecting one or any combination of the following methods to perform any one or several operations on the measured apparent resistivity according to the received processing instructions: single measuring point correction, equal frequency filtering and logarithmic space filtering. Similar to the preprocessing instructions at the original potential difference, the preprocessing instructions here can be the apparent resistivity data of the point number, line number or combination thereof to be processed input by the graphical user interface GUI defined by the computer software program, and / or the preprocessing method selection input after the apparent resistivity anomaly is identified by manual map reading.

[0018] S2 geosteering-based well-seismic combined formation thickness and formation apparent resistivity modeling includes:

[0019] S21, obtaining geological profile data, wherein the geological profile is drawn in combination with existing seismic profiles and geological knowledge;

[0020] S22, performing formation thickness modeling according to the geological profile data and formation apparent resistivity logging data;

[0021] S23, according to the logging elevation data and the apparent resistivity data, the weighted average apparent resistivity, the maximum apparent resistivity and the minimum apparent resistivity corresponding to the formation thickness model are statistically calculated; and the parameters are passed to the formation thickness model: the weighted average apparent resistivity is used as the initial value of the apparent resistivity of the formation thickness model, and the maximum apparent resistivity and the minimum apparent resistivity are used as the upper and lower limits of the model respectively.

[0022] S3's one-dimensional band-source apparent resistivity inversion based on a nonlinear least squares fitting algorithm includes:

[0023] S31, the electric field between the measured electrodes is obtained by forward modeling the formation thickness and apparent resistivity parameters;

[0024] S32, obtain the forward modeled wide-area apparent resistivity by inverse interpolation method based on Emn;

[0025] S33, calculating the residual sum of squares between the forward modeled wide-area apparent resistivity and the measured wide-area apparent resistivity by a least squares fitting function;

[0026] S34, forward iteration, outputs the formation thickness and apparent resistivity that meet the residual square sum requirement as the inversion result.

[0027] As a storage carrier for implementing the present invention, the present invention further provides a computer-readable storage medium on which a computer program is stored, and the computer program implements the steps of the above method when executed by a processor.

[0028] A computer device as a computing carrier for implementing the present invention includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of the above method are implemented.

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

[0030] Seismic data can provide information such as wave impedance interfaces and structures of different layers, but it is difficult to image blank reflection bands with low signal-to-noise ratio in steep and complex areas. Although logging (drilling) data accurately reveals the apparent resistivity information of the formation, it can only evaluate the electrical properties of the drilling position, which is a "one-hole view". Based on geological guidance, the present invention fully integrates and gives full play to the advantages of the three geophysical exploration methods of seismic, logging and electromagnetic, reduces the multi-solution of electromagnetic inversion, improves the inversion depth and accuracy, and provides a strong guarantee for subsequent operations such as profile interpretation, well location deployment, and fluid identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 The flowchart of the electromagnetic inversion method based on well-electric-seismic joint modeling is shown;

[0032] Figure 2 A flow chart of one-dimensional band-source apparent resistivity inversion using a nonlinear least squares fitting algorithm is shown;

[0033] Figure 3 Shows EE x Layout diagram of the device's field source AB and receiving station MN;

[0034] Figure 4 The schematic diagram of N-layer horizontal layered medium in the forward modeling process is shown;

[0035] Figure 5 A comparison diagram before and after potential difference processing-median filtering is shown;

[0036] Figure 6 A comparison diagram of the potential difference treatment before and after the curve smoothing is shown;

[0037] Figure 7 A comparison chart before and after the potential difference treatment-mean regression is shown;

[0038] Figure 8 A comparison diagram before and after the potential difference treatment-curve shift is shown;

[0039] Fig. 9 A comparison diagram before and after potential difference processing-jump point correction is shown;

[0040] Fig.10 A comparison diagram of the current value before and after correction processing is shown;

[0041] Fig.11The comparison diagram of wide-area apparent resistivity processing-single-point correction before and after is shown;

[0042] Fig.12 The comparison diagram of wide-area apparent resistivity processing-equal frequency filtering before and after is shown;

[0043] Fig.13 The comparison diagram of wide-area apparent resistivity processing-logarithmic spatial filtering before and after is shown;

[0044] Fig.14 The inverted initial model stratigraphic thickness map is shown;

[0045] Fig.15 The inverted initial model formation apparent resistivity map is shown;

[0046] Fig.16 The single-point constrained inversion result diagram of formation apparent resistivity is shown;

[0047] Fig.17 A comparison diagram of single-point constrained inversion results of formation apparent resistivity is shown;

[0048] Fig.18 A comparison chart of the single-point inversion value and measured value of wide-area apparent resistivity is shown;

[0049] Fig.19 The section of the inversion result of well-seismic joint modeling is shown. DETAILED DESCRIPTION

[0050] To make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be described clearly and completely in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them.

[0051] The electromagnetic inversion method based on joint well-electric-seismic modeling includes:

[0052] S1, obtain measured apparent resistivity data;

[0053] S2, well-seismic combined formation thickness and formation apparent resistivity modeling based on geosteering;

[0054] S3, one-dimensional band-source apparent resistivity inversion based on nonlinear least squares fitting algorithm;

[0055] S4, output the formation apparent resistivity inversion results.

[0056] This electromagnetic inversion method can be applied to both traditional electromagnetic exploration methods and wide-area electromagnetic exploration methods. Compared with other electromagnetic exploration methods, the difference of the wide-area electromagnetic exploration method is that the calculation method of apparent resistivity is different: the wide-area electromagnetic method defines an exact solution for wide-area apparent resistivity, which is different from the approximate solution of the Kania apparent resistivity of the traditional electromagnetic exploration method. At the same time, since the wide-area electromagnetic method has a greater effective exploration depth than the traditional electromagnetic method, this embodiment uses the wide-area electromagnetic method exploration as an example to illustrate the electromagnetic inversion method based on well-seismic joint modeling.

[0057] refer to Figure 1 Taking the wide-area electromagnetic method as an example, it describes a flow chart of the method from field measured data collection to the final formation apparent resistivity inversion. The field data preprocessing is the process of collecting the original potential difference data from the wide-area receiver, editing, converting and calculating the wide-area apparent resistivity data, which mainly includes the original potential difference data processing, the original current data correction processing and the wide-area apparent resistivity calculation.

[0058] Raw potential difference data processing

[0059] The original potential difference data is the first-hand data collected by the wide-area receiver, which measures the potential difference between two adjacent electrodes MN. In the actual field data collection process, affected by various factors such as high-voltage lines, large towns, and wiring firmness, the original potential difference data will inevitably have multiple interferences such as random interference and fixed interference, which are mainly manifested in the distortion of the potential difference corresponding to some frequency points of the measuring point and the up and down deviation of the potential difference curve corresponding to the overall frequency point. If the original potential difference data is used directly for wide-area apparent resistivity calculation, unnecessary deviations or even errors will inevitably occur and affect the reliability of the final formation apparent resistivity inversion. Therefore, the main preprocessing methods for the quality deviation of the original potential difference data caused by multiple factors include one or any combination of the following methods: median filtering, mean filtering, mean regression, curve translation, and jump point correction.

[0060] Each processing method is suitable for different processing scenarios. Any combination of processing methods can meet the impact of multiple random / fixed interference factors on the original potential difference data. The detailed introduction of various pre-processing methods is as follows:

[0061] (1) Median filtering: It is mainly used to eliminate the distortion points caused by the abnormal increase of the potential difference value of individual frequency points due to interference such as high-voltage lines, and can correct them to ensure the continuity of the potential difference curve. Figure 5The figure shows the comparison before and after median filtering. Before median filtering, there are two jump points in the high frequency part (in the blue dotted box in the left figure). The jump points are caused by the presence of high-voltage transmission lines above 110KV or large towns, which cause abnormal increase in the potential difference of individual frequency points (this type of interference will increase the potential difference, not decrease it). After median filtering, the jump points in the high frequency part disappear, and there is no obvious impact on the surrounding data (in the blue dotted box in the right figure). Median filtering can effectively remove noise and retain the edge information of the signal, and is suitable for denoising obvious distortion points.

[0062] (2) Curve smoothing, also known as mean filtering: It is used to eliminate the phenomenon that the curve in the potential difference data is not smooth and continuous. It can ensure the continuity of the potential difference curve and make it smoother. Compared with the median filter, it has advantages in processing single distortion points with large amplitude changes. Curve smoothing is more targeted at processing a local weak discontinuity. Figure 6 Comparison before and after the curve smoothing process. There is a weak discontinuity in the potential difference in the intermediate frequency part (in the blue box on the left), which is generally caused by random interference caused by environmental factors (such as wind, etc.), which can be eliminated by curve smoothing. After processing, the curve is smooth and continuous, and the morphological characteristics are more obvious (in the blue box on the right, the smoothing method is 5-point mean smoothing). Mean filtering can reduce the impact of noise, but it will blur the edge of the signal, which is suitable for denoising random slight weak interference.

[0063] (3) Mean regression: If the overall data quality of the current measuring point is poor, but the data quality of the previous and next measuring points is good, the shape is basically the same, and the amplitude difference is very small, then for the processing of the current measuring point, the average value of the potential difference of the previous and next measuring points can be used to replace the potential difference value of the current measuring point, which is equivalent to taking a two-point average of the current measuring point (excluding the current measuring point). Figure 7 Comparison before and after processing to update the average value: before processing (left picture), the overall data quality is poor (red), but the shape is basically consistent with the previous and next measuring points (blue and green) and the amplitude is close. Due to the lateral continuity of the formation, the average value of the previous and next two measuring points is used to replace the current measuring point (right picture).

[0064] (4) Curve translation (upward and downward): Usually, the potential difference curve is relatively continuous and does not have large step jumps. If such jumps occur, it is necessary to move the step jump segment of the current measuring point data up or down according to the data of the previous and next measuring points to return the step jump to its original position and ensure the continuity of the curve. Figure 8 Comparison before and after the downward shift: The potential difference amplitude of the low-frequency part increases abnormally as a whole, which is discontinuous with the previous data, and a step-like jump occurs (the blue rectangular area on the left). By moving the low-frequency part downward as a whole, the step-like jump is eliminated and the data continuity is ensured (the blue rectangular area on the right).

[0065] (5) Jump point correction: When there are many distortion points in the potential difference data and other processing methods are not applicable, the potential difference values ​​at different frequency points of the measuring point can be corrected by single / multi-jump point correction. The purpose is also to ensure that the processed curve retains the original curve shape as a whole and is smooth and continuous. Fig. 9 This is a comparison chart before and after jump point correction: the high-frequency and low-frequency parts with poor data quality were processed, and the mid-frequency part was retained (right picture) to ensure that the curve shape is smooth after processing while respecting the original data as much as possible.

[0066] Raw current data correction processing

[0067] Wide-area electromagnetic generators generally transmit current in the order of high-medium-low frequency groups. In order to ensure the stability of the transmitted current, the transmission voltages of different frequency groups are different: for the high-frequency group, as the frequency decreases, the voltage needs to be continuously increased; for the medium-frequency group, the voltage required to transmit the same current is extremely high, and it needs to be guaranteed to be above 1000V; the same is true for the medium-frequency-low frequency connection. Therefore, it is difficult for the current at the connection between the high-medium frequency group and the medium-low frequency group to transition smoothly, and jumps will occur. However, the current data must be used as an important parameter for the subsequent wide-area apparent resistivity calculation. The current at the frequency group connection must be corrected to ensure that the current curve is smooth and continuous, thereby ensuring that the results of the wide-area apparent resistivity calculation are continuous and reliable, without discontinuous distortion points.

[0068] By selecting the frequency point, select the connection of different frequency groups (high-medium or medium-low) or other current jump point positions for smoothing to ensure smooth and continuous frequency group connection and the entire current curve. Typical current correction methods such as five-point mean smoothing actually require only editing and processing the current curve of a single measuring point when smoothing the current data of the measuring line, and then extending it to the remaining measuring points, that is, all measuring points of the same source use the same current data. Fig.10 The figure shows the comparison of current data before and after correction. The left and right figures are before and after current correction, respectively. It can be seen that the current continuity at the connection of high-frequency group-medium-frequency group and medium-frequency group-low-frequency group (shown in two green circles) has been improved, and the local weak discontinuity (medium-frequency group and low-frequency group) has also been smoothed to ensure the reliability of subsequent wide-area apparent resistivity calculation; and the same current data is used for the same survey line, which also ensures the uniformity of subsequent wide-area apparent resistivity calculation.

[0069] Wide-area apparent resistivity calculation

[0070] Before calculating the wide-area apparent resistivity, it is necessary to sort the coordinates first and merge the field source coordinates, measuring point coordinates, measuring point frequency, measuring point current and measuring point potential difference data. Usually, they are spliced ​​and merged in the order of point number, field source A pole x coordinate (Ax), field source A pole y coordinate (Ay), field source B pole x coordinate (Bx), field source B pole y coordinate (By), receiving electrode M pole x coordinate (Mx), receiving electrode M pole y coordinate (My), receiving electrode N pole x coordinate (Nx), receiving electrode N pole y coordinate (Ny), frequency, current, potential difference and error.

[0071] After the coordinates are sorted, the wide-area apparent resistivity is calculated by the inverse interpolation method (also called the inverse iteration method). The formula principle is as follows:

[0072] mf=1,2,…,pd

[0073]

[0074] rho_range = logspace(0,5,51)

[0075] rho_WFEM(mf)=interp1(Emn_p(mf),rho_range,Emn(mf))

[0076] Where rho_WFEM represents the wide-area apparent resistivity to be calculated, interp1 represents one-dimensional interpolation (generally cubic spline interpolation is used), rho_range is the given wide-area apparent resistivity sequence, Emn_p is the electric field between the MN electrodes calculated based on rho_range, and Emn is the measured electric field between the MN electrodes, which is calculated based on the measured potential difference Vmn ( is the distance between the receiving electrodes MN), mf represents the frequency point sequence, and pd represents the number of frequency points.

[0077] logspace represents a logarithmic interval generator, which generates 51 points between 0 and 5 with a base of 10 (interval 0.1), that is, the range of the general wide-area apparent resistivity series is 10 0 Up to 10 5 (10,000 to 100,000). It can be seen that when other variable parameters are determined, only Emn_p needs to be calculated to interpolate the required wide-area apparent resistivity. x The device requires that the source AB is parallel to the receiving station MN, and only a single E x Component, and obtain Emn_p by taking the absolute value.

[0078] i=1,2,...,51

[0079]

[0080] Emn_p(i)=abs(E x )

[0081] In the above formula, PE represents the electric dipole, I represents the corrected current data, is the distance between the field sources AB, pi is a constant, and the value of pi is 3.141592653. r represents the distance between the midpoint of the field source AB and the midpoint of MN, It can be understood as the angle between the direction of the field source AB and the line connecting the midpoint of AB to the midpoint of MN.

[0082] i represents an imaginary number, k represents a wave number, and the expression is as follows:

[0083]

[0084] w is the angular frequency, and w = 2*pi*f, f is the frequency (unit Hz), μ is the air magnetic permeability, the value is 4*pi*10 -7 ; ε represents the dielectric constant of vacuum, which is 8.854*pi*10 -12 In actual use, the wave number calculation formula can be simplified, because the former is much smaller than the latter and its influence on the result is negligible. The simplified wave number k calculation formula commonly used is as follows:

[0085]

[0086] After calculating the measured apparent resistivity, apparent resistivity processing is also performed, which includes any one or more of the following operations: single measuring point correction, equal frequency filtering and logarithmic space filtering.

[0087] (1) Single point calibration

[0088] The wide-area apparent resistivity data is calculated using the pre-processed potential difference data, so it is relatively smooth and continuous, and only requires local fine-tuning. Over-processing will affect the authenticity of subsequent apparent resistivity inversion. Single-point correction can eliminate the influence of static displacement to a certain extent. Fig.11 As shown in the figure, the calculated wide-area apparent resistivity is smooth and interference-free as a whole. It should be noted that the apparent resistivity of the high-frequency first branch or the low-frequency tail branch of some measuring points will show a step-like jump (sudden change and increase upward, as shown in the blue rectangular box on the left). The abnormal part can be moved downward as a whole to ensure the gradual change and continuity of the apparent resistivity of the first and tail branches.

[0089] (2) Equal frequency filtering

[0090] Unlike single-point correction, which processes by measuring point, equal-frequency filtering (also called equal-frequency curve editing) performs lateral processing based on frequency points. The purpose is to improve the lateral continuity of the subsequent inversion profile, and to eliminate the influence of static displacement to a large extent. Considering the relationship between independent variables and dependent variables, in layman's terms, single-point correction uses a single measuring point as the independent variable, and the wide-area apparent resistivity of all frequency points of this measuring point is the dependent variable; while equal-frequency filtering uses the frequency point as the independent variable, and the apparent resistivity of all measuring points of this frequency point is the dependent variable. Equal-frequency filtering can achieve targeted and fine processing by box-selecting the measuring section, or it can achieve one-click processing of all frequency points by batch smoothing (usually using 5-point mean smoothing). Fig.12 This is a comparison before and after equal-frequency filtering (taking 256 Hz single frequency point as an example). The horizontal axis is the point number, and the vertical axis is the logarithm of the wide-area apparent resistivity with a base of 10. The filtering effect is shown in the two blue dotted rectangular boxes. Before filtering, there are more "roughness" in the horizontal direction, and the amplitudes of many measuring points are discontinuous with the previous and next measuring points. After filtering, the horizontal direction is smooth and continuous, eliminating some static effects.

[0091] (3) Logarithmic spatial filtering

[0092] Single-point correction processes the measurement point as the smallest unit, while equal-frequency filtering processes the frequency point as the smallest unit, both of which are one-dimensional filtering. Logarithmic space filtering is a two-dimensional filtering process for the apparent resistivity file, aiming to improve the spatial (horizontal and vertical) continuity of the apparent resistivity profile. Logarithmic space filtering is a two-dimensional space (2D) filtering, that is, it considers both horizontal and vertical directions at the same time, the horizontal axis is the point number, the vertical axis is the logarithmic frequency with a base of 10, and the color scale shows the logarithmic wide-area apparent resistivity with a base of 10. Fig.13 The "frequency-apparent resistivity" pseudo-sections are original, equal frequency filtering, and logarithmic space filtering. The horizontal axis is the point number, the vertical axis is the logarithmic frequency, and the color scale represents the logarithmic wide-area apparent resistivity (all with a base of 10). It can be seen that the original wide-area apparent resistivity single-point curve is relatively smooth and continuous, but its spatial distribution characteristics are relatively rough (in the red rectangular frame), and the horizontal and vertical continuity is very poor; after equal frequency filtering, the spatial distribution is significantly improved, the horizontal and vertical (especially the horizontal) continuity is significantly enhanced, the electrical characteristics are obvious, but there is still a slight discontinuity in the vertical direction; after logarithmic space filtering, the spatial characteristics are further clear, and the slight discontinuity in the vertical direction is improved.

[0093] S2 geosteering-based well-seismic combined formation thickness and formation apparent resistivity modeling includes:

[0094] S21, draw geological profiles based on existing seismic profiles and geological knowledge;

[0095] S22, modeling the formation thickness based on geological profiles and formation apparent resistivity logging data;

[0096] S23, according to the logging elevation data and the apparent resistivity data, the weighted average apparent resistivity, the maximum apparent resistivity and the minimum apparent resistivity corresponding to the formation thickness model are statistically calculated; the parameters are passed to the formation thickness model, the weighted average apparent resistivity is used as the initial value of the apparent resistivity of the formation thickness model, and the maximum apparent resistivity and the minimum apparent resistivity are used as the upper and lower limits of the model respectively.

[0097] Seismic data provides relatively detailed stratigraphic and structural information, which can provide a reference for geological profile drawing. However, for areas such as high-steep structural belts where the quality of seismic data is poor, stratigraphic inference is carried out in combination with surface geological outcrops and related geological knowledge, and geological profiles are drawn, making full use of existing seismic data and geological knowledge as the basis for stratigraphic thickness modeling. The length of the drawn geological profile is consistent with the actual survey line length of the wide-area electromagnetic method, the surface should be consistent with the actual measured point coordinates, the horizontal coordinate is the point number (mileage), the vertical coordinate is the elevation (unit: m), and the profile should contain information such as stratigraphic boundaries and stratigraphic codes.

[0098] According to the geological profile and logging apparent resistivity data, the formation thickness modeling is carried out by means of elevation setting (editing). The first layer elevation is the surface elevation (maximum elevation) of all measuring points, the last layer is the minimum elevation of inversion, and the elevation of the middle large layer is set (edited) according to the formation boundary line. For thicker large layers or large layers with obvious changes in apparent resistivity on the geological profile, in the preferred implementation method, multiple small layers are finely divided within the large layer according to the logging data. After the elevation setting (editing) is completed, the first-order derivative of the elevation of each measuring point is obtained to obtain all the formation thickness model information of the measuring point. In practice, the typical number of inversion layers is 10 to 30 layers. Fig.14 It is a stratigraphic thickness map of the inversion initial model jointly established based on geological steering and well seismic. The horizontal axis is the point number, the vertical axis is the elevation (unit: meter), the first elevation line from top to bottom is the surface elevation of the measuring point, and the last elevation is the minimum inversion elevation (-6500 meters). The intermediate elevation lines are based on the wave impedance interface of seismic data, supplemented by geological knowledge and logging apparent resistivity.

[0099] The geological model consists of two parts: formation thickness and formation apparent resistivity. After the formation thickness model is established, the formation apparent resistivity attributes are filled in to perform formation apparent resistivity modeling. First, the weighted average apparent resistivity, minimum apparent resistivity and maximum apparent resistivity corresponding to the formation thickness model are statistically calculated based on the two columns of well logging data elevation and apparent resistivity, and their parameters are passed to the geological model. The weighted average apparent resistivity is used as the initial value of the geological model apparent resistivity, and the statistical maximum apparent resistivity and minimum apparent resistivity are used as the upper and lower limits of the model, respectively. Among them, the weighted average apparent resistivity is the ratio of the sum of the product of the thickness of each well logging layer and its apparent resistivity value in the single layer of the thickness model to the thickness of the single layer. Fig.15It represents the apparent resistivity map of the inversion initial model formation established based on the well logging apparent resistivity calibration. The inversion initial model formation apparent resistivity is the weighted average apparent resistivity calculated based on the well logging apparent resistivity curve and the inversion initial model layer thickness.

[0100] S3's one-dimensional band-source apparent resistivity inversion based on a nonlinear least squares fitting algorithm includes:

[0101] S31, the electric field Emn between the measured MN electrodes is obtained by forward modeling the formation thickness and apparent resistivity parameters;

[0102] S32, according to Emn, the forward modeled wide-area apparent resistivity rho_mn(mf) is obtained by the inverse interpolation method;

[0103] S33, calculating the residual sum of squares between the forward modeled wide-area apparent resistivity rho_mn(mf) and the measured wide-area apparent resistivity rho_WFEM(mf) by a least squares fitting function;

[0104] S34, forward iteration, outputs the formation thickness and apparent resistivity that meet the residual square sum requirement as the inversion result.

[0105] Forward modeling is the basis of inversion. Inversion is completed through uninterrupted forward modeling iterations. It is the process of obtaining the best geological model by back-calculating the measured wide-area apparent resistivity. The process of calculating the measured wide-area apparent resistivity rho_WFEM(mf) using the electric field electrode difference Emn has been introduced in detail in the previous article. The following only introduces the process of forward modeling the electric field from parameters such as formation apparent resistivity and formation thickness.

[0106] The electric dipole source (field source A and B) is to supply harmonic current to the grounded conductor of finite length, and to establish an electromagnetic field of corresponding frequency underground by changing the frequency of the supply current. The resistivity and thickness of the nth layer in the N-layer horizontal layered medium are denoted as ρn and hn, respectively. Figure 3 .

[0107] The center of the electric dipole is taken as the coordinate origin, the x-axis points to the direction of the dipole moment, and the z-axis points vertically downward. Using the conversion relationship between cylindrical coordinates and rectangular coordinate systems, the expression of the surface electromagnetic field component in the rectangular coordinate system under quasi-static conditions is obtained:

[0108]

[0109] In the formula, R and R * It is a function that connects the electrical parameters of each layer of the formation; IdL is the dipole moment (the product of the current and the length of the field source AB); is the angle between the dipole moment direction and the line from the dipole center to the measuring point; r is the transmitting and receiving distance (the distance from the midpoint of the field source AB to the midpoint of MN); w is the angular frequency, and w = 2πf, μ is the magnetic permeability in the air, and μ = 4*pi*10-7 ; The physical meaning of m is a space constant, which has the dimension of the inverse of distance; k i is the wave number of the i-th electrical layer; ρ i is the resistivity of the i-th electrical layer; h i is the thickness of the i-th electrical layer; J0(mr) and J1(mr) are the zero-order Bessel function and the first-order Bessel function with mr as the variable, respectively;

[0110] For the E-Ex device, first forward model to obtain E x Component, let Emn = abs(E x ), and then the inverse interpolation method is used to obtain the final forward modeled wide-area apparent resistivity.

[0111] If the measured wide-area apparent resistivity is rho_WFEM and the forward-modeled wide-area apparent resistivity is rho_mn, the least squares fitting inversion can be achieved by establishing the following function:

[0112]

[0113] pd represents the number of frequency points, log represents the logarithm with base e, and FF is the function value. If the measured wide-area apparent resistivity rho_WFEM(mf) of each frequency point is basically consistent with the forward modeled wide-area apparent resistivity rho_mn(mf), their ratio is close to 1, and the function value FF(mf) should be infinitely close to 0. The formation thickness and apparent resistivity at this time are output as the inversion results.

[0114] Changes in the two attributes of formation thickness and apparent resistivity will affect the calculation of forward wide-area apparent resistivity. If both are considered at the same time, the inversion process will be extremely cumbersome, the number of inversion iterations will be huge, and the amount of calculation will increase significantly. Since the layer thickness attribute has been determined based on geological, seismic and logging data during the joint well-seismic modeling, in the actual inversion process, it is preferred to only use the formation apparent resistivity attribute as the variable to be inverted, and obtain the best formation apparent resistivity inversion result through iterative optimization.

[0115] The time consumption table of single-point constrained inversion is as follows. Fast forward optimization iteration is the key to improving inversion efficiency and accuracy.

[0116]

[0117] Forward optimization methods include but are not limited to the following aspects:

[0118] ① In forward modeling, the usual practice is to set PE (electric dipole parameter) to 1. Here we will strictly follow Formula, calculation is performed based on the measured current to improve the accuracy of the result;

[0119] ②During forward modeling, it is necessary to calculate R and R* Parameters, for multilayer models R and R * The electrical parameters of each layer are deduced from bottom to top. The usual approach is to use the frequency point, number of layers, and filter coefficients (i.e. R and R * The electrical parameter iteration coefficients of the adjacent lower layer to the upper layer in the expression of ) are calculated by triple loops. The general process is: frequency point loop; for each frequency point, the forward modeling layer loop; for each forward modeling layer, the filter coefficient loop. R and R of this type of triple loop * The parameter solution method will undoubtedly greatly increase the amount of calculation. As a preferred implementation method, according to R and R * The expression of the parameters is calculated in the form of a two-dimensional matrix: the calculation of R and R* includes the forward frequency loop. For each frequency point, the parameters R and R* to be determined are calculated in the form of a two-dimensional matrix. The rows and columns of the matrix represent the filter coefficient sequence and the forward layer sequence, respectively. The R and R* of the bottom layer are set to 1. Starting from the second to last layer, the R and R* functions of each layer of each frequency point are calculated from bottom to top by cross-row point multiplication and iteration. In this way, the triple loop is reduced to one dimension, reducing the loop level of the algorithm program implementation process to improve the execution efficiency of the program code, thereby improving the final inversion speed;

[0120] ③ The Ex expression for forward modeling of apparent resistivity contains zero-order and first-order Bessel functions. Since the Bessel function is not only a high-frequency oscillation function, but also a slowly decaying function, these characteristics will affect the accuracy of the calculation. The Hankel integral can only obtain an analytical solution under uniform half-space conditions, but generally no analytical solution can be obtained in other cases, so only numerical integration methods or digital filtering methods can be used. As a preferred embodiment, the Hankel integral is calculated using 61-point Hankel J0 transform linear filter coefficients and 61-point J1 transform linear filter coefficients. While ensuring the calculation speed, better forward modeling results can be obtained compared to other filter coefficients.

[0121] Before optimization, the average time consumption of single-point constrained inversion was 40.561s, and after optimization, the average time consumption was 10.985s, and the calculation efficiency was improved by about 400%. Considering that the current wide-area electromagnetic survey line is generally within 20km in length, if the survey line length is 20km, the point spacing is 100m, and the single-point inversion time is 11s, the number of survey points is 201, and the total inversion time does not exceed 37min. If the conventional 4-core computer is used for parallel inversion, the total time consumption will be controlled within 10min, which is completely acceptable for the well-seismic joint inversion proposed in the present invention.

[0122] Taking the typical MATLAB toolbox as an example, the nonlinear least squares optimization solution is introduced. The MATLAB toolbox provides the lsqnonlin function to solve the nonlinear least squares optimization problem. Its calling format is X=lsqnonlin(fun,x0,lb,ub,options).

[0123] X is the apparent resistivity of the formation obtained by the final inversion, which is consistent with the initial layer thickness or the number of apparent resistivities; fun represents the least squares fitting function, which consists of three parts:

[0124] ①Emn is obtained by forward modeling parameters such as formation thickness and apparent resistivity;

[0125] ②Use Emn and inverse interpolation method to calculate the forward wide-area apparent resistivity;

[0126] ③ Calculate the residual sum of squares of the forward modeling and measured wide-area apparent resistivity through the FF function;

[0127] x0 is the initial value of apparent resistivity of the established geological model, lb and ub are the lower and upper limits (minimum and maximum values) of apparent resistivity of the formation model respectively;

[0128] options is the optimization option. This function provides two algorithms for optimization calculation: trust-region-reflective algorithm and levenberg-marquardt algorithm. Here we use the trust-region-reflective algorithm because the levenberg-marquardt algorithm is an unconstrained algorithm. It only needs to specify the initial value x0 and does not require lower and upper limits. The trust-region-reflective algorithm is a constrained algorithm. Its inversion starts from the initial value x0, and the final inversion result will be taken within the given lb and ub.

[0129] In addition, during inversion, multi-core parallel inversion is started using the parfor loop to invert multiple measurement points at the same time, greatly shortening the inversion time.

[0130] After the inversion is completed, the inversion model is converted into three columns: point number, elevation and formation apparent resistivity, to facilitate subsequent detailed mapping in data visualization tools (such as Surfer software).

[0131] Fig.16The figure shows the result of single-point constrained inversion of formation apparent resistivity. The horizontal axis is the apparent depth (unit: meter), and the vertical axis is the inverted apparent resistivity (base 10 logarithm, unit: ohm meter). According to the constraints of geosteering and well-seismic joint modeling, the number of inversion layers is determined to be 15, with a total of 4 curves. The green and blue lines represent the given formation apparent resistivity thresholds (inversion minimum and maximum values), the pink solid line represents the initial value of the formation apparent resistivity (weighted average apparent resistivity of different formations in the logging curve), and the black line represents the final inversion value of the formation apparent resistivity. It can be seen that after the initial value and upper and lower limits are given, the nonlinear least squares fitting algorithm will automatically optimize within the given range to achieve the best inversion effect.

[0132] Attached Fig.17 and attached Fig.18 This is a single point apparent resistivity inversion comparison.

[0133] Attached Fig.17 The black line represents the constrained inversion result (i.e. Fig.16 The red line represents the final result of the inversion in ); the red line represents the unconstrained inversion result. For comparison, the unconstrained inversion also sets the strata to 15 layers (of equal thickness). From the perspective of electrical characteristics, the unconstrained inversion only inverts four changes: low-high-low-high, while the constrained inversion inverts a total of 10 electrical characteristics: low-high-low-high-low-high-low-high-low-high-low-high. Overall, both curves show the characteristics of low-high-low-high variation, but because the constrained inversion integrates geological, seismic, and logging information at the same time, the inversion results are more refined and more in line with geological laws.

[0134] Attached Fig.18 This is a comparison chart of the single-point inversion and measured values ​​of wide-area apparent resistivity.

[0135] The black line is the measured wide-area apparent resistivity, the blue line is the unconstrained inversion wide-area apparent resistivity, and the red line is the constrained inversion wide-area apparent resistivity. In essence, the wide-area apparent resistivity is the forward modeling result of the formation apparent resistivity. Although the final inversion curve is inevitably different from the measured curve, it is not difficult to see that the constrained inversion wide-area apparent resistivity is closer to the measured value than the unconstrained one (as shown in the black dotted rectangle).

[0136] Attached Fig.19 The final result profile is inverted for well-seismic joint modeling.

[0137] It can be seen that compared with Fig.14 and 15 The initial inversion model is attached. Fig.19The final result of the inversion effectively reduces the multi-solution of the inversion due to the combination of seismic geology, electromagnetic and logging information. The lateral variation characteristics of the profile are obvious. For the two sets of low-resistance layers in the profile, the apparent resistivity is relatively high in the high-steep structural belt (in the red dotted box), and relatively low in the low-rise (shown in the blue dotted box); vertically, a total of 8 main electrical layers are divided, including low-high-low-high-low-high-low-high-low-high, and the stratification characteristics are obvious. The apparent resistivity distribution of the profile is consistent with the geological law in both the horizontal and vertical directions, and the minimum inversion elevation can reach -5000 meters or even lower, which is difficult to achieve with unconstrained inversion.

[0138] Therefore, the following detailed description of the embodiments of the present invention is not intended to limit the scope of the invention claimed for protection, but merely represents some 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.

[0139] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features and technical solutions in the embodiments may be combined with each other.

[0140] The above embodiments are only used to illustrate the present invention and are not intended to limit the technical solutions described in the present invention. Although the present invention has been described in detail in this specification with reference to the above embodiments, the present invention is not limited to the above specific implementation methods. Any technical solutions and improvements that modify or replace the present invention without departing from the spirit and scope of the invention are all included in the scope of the claims of the present invention.

Claims

1. The electromagnetic inversion method based on well-electric-seismic joint modeling is characterized by: include: S1, obtain measured apparent resistivity data; S2, well-seismic combined formation thickness and formation apparent resistivity modeling based on geosteering; S3, one-dimensional band-source apparent resistivity inversion based on nonlinear least squares fitting algorithm; S4, output the formation apparent resistivity inversion results.

2. The method according to claim 1, characterized in that S2 geosteering-based well-seismic combined formation thickness and formation apparent resistivity modeling includes: S21, obtaining geological profile data, wherein the geological profile is drawn in combination with existing seismic profiles and geological knowledge; S22, performing formation thickness modeling according to the geological profile data and formation apparent resistivity logging data; S23, according to the logging elevation data and the apparent resistivity data, the weighted average apparent resistivity, the maximum apparent resistivity and the minimum apparent resistivity corresponding to the formation thickness model are statistically calculated; and the parameters are passed to the formation thickness model: the weighted average apparent resistivity is used as the initial value of the apparent resistivity of the formation thickness model, and the maximum apparent resistivity and the minimum apparent resistivity are used as the upper and lower limits of the model respectively.

3. The method according to claim 1, characterized in that The one-dimensional band-source apparent resistivity inversion based on the nonlinear least squares fitting algorithm of S3 only takes the formation apparent resistivity as the variable to be inverted, and the formation thickness is determined by the well-seismic combined formation thickness model based on geosteering.

4. The method according to any one of claims 1 to 3, characterized in that: The measured apparent resistivity data obtained by S1 is the apparent resistivity data calculated by wide-area electromagnetic method, further including: S11, performing potential difference preprocessing on the original potential difference data collected by the wide-area electromagnetic method receiver; S12, preprocessing the current value of the transmission current at the non-smooth current position of the wide-area electromagnetic method transmitter; S13, the measured apparent resistivity value is calculated using the inverse interpolation method.

5. The method according to claim 4, characterized in that The potential difference preprocessing described in S11 includes correcting the original potential difference abnormal points or areas according to the received preprocessing instruction by selecting one of the following methods or any combination of several methods: median filtering, mean filtering, mean regression, curve translation and jump point correction.

6. The method according to claim 4, characterized in that The current value preprocessing described in S12 includes smoothing the connection between different frequency groups of a single measuring point or other non-smooth emission current value abnormal points or areas according to the received preprocessing instructions, so as to make the frequency group connection and the complete current curve smooth and continuous; and, extending the complete current curve of the single measuring point to the remaining measuring points of the same field source.

7. The method according to claim 4, characterized in that Before calculating the measured apparent resistivity, S13 first sorts the coordinates and then uses the following formula to calculate the measured apparent resistivity rho_WFEM: rho_WFEM(mf)=interp1(Emn_p(mf),rho_range,Emn(mf)) rho_range = logspace(0,5,51) mf=1,2,…,pd Where interp1 represents one-dimensional interpolation, rho_range is a given wide-area apparent resistivity series, Emn_p is the electric field between MN electrodes calculated based on rho_range, logspace represents a logarithmic interval generator that generates 51 points between 0 and 5 with a base of 10, Emn is the measured electric field between MN electrodes, and Vmn is the measured potential difference. is the distance between the receiving electrodes MN, mf is the frequency point sequence, and pd is the number of frequency points.

8. The method according to claim 7, characterized in that The original potential difference data collected by the wide-field electromagnetic receiver is collected in a way that the receiving electrode MN is parallel to the field source AB. Emn_p is obtained using the following formula: Emn_p(i)=abs(E x ) i＝1,2,…,51 Where PE represents the electric dipole, I represents the corrected current data, is the distance between the field sources AB, pi is the pi constant, r is the distance between the midpoint of the field source AB and the midpoint of MN, It represents the angle between the direction of the field source AB and the line connecting the midpoint of AB to the midpoint of MN; i represents an imaginary number, k represents a wave number, and the expression is as follows: Among them, w=2*pi*f, represents angular frequency, f represents frequency, μ represents air magnetic permeability, and ε represents vacuum dielectric constant.

9. The method according to claim 8, characterized in that After calculating the measured apparent resistivity, S13 also performs apparent resistivity processing, which includes selecting one or any combination of the following methods to perform any one or several operations on the measured apparent resistivity according to the received processing instructions: single measuring point correction, equal frequency filtering and logarithmic space filtering.

10. The method according to claim 9, characterized in that S3's one-dimensional band-source apparent resistivity inversion based on a nonlinear least squares fitting algorithm includes: S31, the electric field Emn between the measured MN electrodes is obtained by forward modeling the formation thickness and apparent resistivity parameters; S32, according to Emn, the forward modeled wide-area apparent resistivity rho_mn(mf) is obtained by the inverse interpolation method; S33, calculating the residual sum of squares between the forward modeled wide-area apparent resistivity rho_mn(mf) and the measured wide-area apparent resistivity rho_WFEM(mf) by a least squares fitting function; S34, forward iteration, outputs the formation thickness and apparent resistivity that meet the residual square sum requirement as the inversion result.

11. The method according to claim 10, characterized in that S31 The electric field Emn between the measured MN electrodes obtained by forward modeling the formation thickness and apparent resistivity parameters includes: In the formula, R and R * It is a function that connects the electrical parameters of each layer of the formation; IdL is the dipole moment (the product of the current and the length of the field source AB); is the angle between the dipole moment direction and the line from the dipole center to the measuring point; r is the transmitting and receiving distance (the distance from the midpoint of the field source AB to the midpoint of MN); w is the angular frequency, and w = 2πf, μ is the magnetic permeability in the air, and μ = 4*pi*10 -7 ; The physical meaning of m is a space constant, which has the dimension of the inverse of distance; k i is the wave number of the i-th electrical layer; ρ i is the resistivity of the i-th electrical layer; h i is the thickness of the i-th electrical layer; J0(mr) and J1(mr) are the zero-order Bessel function and the first-order Bessel function with mr as the variable, respectively; For the E-Ex device, first forward model to obtain E x Component, let Emn = abs(E x ).

12. The method according to claim 11, characterized in that The calculation of R and R* includes the forward frequency cycle. For each frequency point, the parameters R and R* are calculated in the form of a two-dimensional matrix. The rows and columns of the matrix represent the filter coefficient sequence and the forward layer sequence, respectively. The R and R* of the bottom layer are set to 1. Starting from the second to last layer, the R and R* functions of each frequency point and each layer are calculated from bottom to top by cross-row dot multiplication and iteration.

13. The method according to claim 11, characterized in that The Hankel integral is calculated using the 61-point Hankel J0 transform linear filter coefficients and the 61-point J1 transform linear filter coefficients.

14. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 13 are implemented.

15. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 13 are implemented.

Images