A VSP first arrival travel time tomography method and system

Abstract

The present invention relates to a VSP first arrival travel time tomography method and system, and belongs to the technical field of marine mineral resource exploration and seabed stratum detection. First, the ground seismic profile is processed to obtain ground seismic structure information data to construct an initial VSP velocity model. Then, the first arrival travel time calculation value is calculated according to the initial VSP velocity model, and the first arrival travel time measurement value and the first arrival travel time calculation value are subtracted to obtain a time residual matrix. The time residual matrix is ​​used as input to solve the inversion equation to obtain an updated VSP velocity model until the updated VSP velocity model meets the predetermined accuracy, thereby using the ground seismic structure information data as prior knowledge to constrain the inversion process of the VSP velocity model, and improve the accuracy of the VSP first arrival travel time tomography.

Description

Technical Field

[0001] The present invention relates to the technical field of marine mineral resource exploration and seabed stratum detection, and in particular to a VSP first-arrival traveltime tomography method and system under the constraint of ground seismic structure information data. Background Art

[0002] Seismic exploration refers to the propagation of seismic waves caused by the excitation of artificial seismic sources in underground media. Due to the differences in elasticity and density of underground media, transmission, reflection, refraction and other phenomena occur, carrying the layer and lithology information of the underground media to the surface. Geophysicists process and interpret seismic wave signals, analyze the propagation law of seismic waves in underground media and further infer the properties and morphology of underground rock formations. There are many methods of seismic exploration. According to the distribution of seismic sources and detectors, they can be divided into: ground seismic exploration method, vertical seismic profile (VSP) method, and cross-well seismic method. Ground seismic refers to onshore ground, marine streamer and seabed node earthquakes. The source and the detector of the ground seismic exploration method are approximately in the same horizontal direction (hereinafter referred to as the horizontal plane). The source is arranged on the ground or below the sea surface (streamer), and the detector is arranged on the ground, below the sea surface (streamer) or on the seabed. The VSP method is an exploration method in which seismic waves are excited by a horizontal source and the detector fixed on the well wall receives the seismic records. Compared with surface seismic methods, VSP seismic exploration methods have higher resolution and signal-to-noise ratio, can provide spatial and temporal information of the strata near the wellbore, and describe the seismic structure of the strata near the well in more detail.

[0003] The observation method of VSP seismic data reduces travel time and ray propagation route, reduces energy attenuation caused by near-surface low-velocity belt, has high seismic record quality, is easy to pick up the initial arrival travel time, and tomography technology is an effective way to establish VSP velocity model. However, the amount of VSP seismic data is less than that of ground seismic data, the range of detection points is small, and the target imaging range is limited to the well area. Therefore, the comprehensive use of prior information constraints to improve inversion accuracy is still the key direction of VSP seismic data inversion research, and it is also a difficult problem that needs to be solved urgently in production practice.

[0004] Based on this, there is an urgent need for a VSP first-arrival traveltime tomography method and system that can improve the inversion accuracy. Summary of the invention

[0005] The purpose of the present invention is to provide a VSP first-arrival traveltime tomography method and system, which utilizes ground seismic structure information data as prior knowledge to constrain the inversion process of the VSP velocity model and improve the accuracy of VSP first-arrival traveltime tomography.

[0006] To achieve the above object, the present invention provides the following solutions:

[0007] A VSP first arrival travel time tomography method, the tomography method comprising:

[0008] Processing the ground seismic profile to obtain ground seismic structure information data, and constructing an anisotropic diffusion equation based on the ground seismic structure information data;

[0009] Taking one-dimensional velocity data as input, solving the anisotropic diffusion equation to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data;

[0010] Extract first arrival traveltime measurements from VSP seismic data;

[0011] According to the initial VSP velocity model, the first arrival travel time calculation value is calculated by using the fast scanning method;

[0012] Subtracting the first arrival travel time measurement value from the first arrival travel time calculation value to obtain a time residual matrix;

[0013] Taking the time residual matrix as input, solving the inversion equation to obtain an updated VSP velocity model;

[0014] Determining whether the updated VSP velocity model meets a predetermined accuracy;

[0015] If yes, the updated VSP velocity model is used as the final VSP velocity model;

[0016] If not, the updated VSP velocity model is used as the initial VSP velocity model in the next iteration, and the process returns to the step of "calculating the first arrival travel time by using the fast scanning method according to the initial VSP velocity model".

[0017] A VSP first arrival travel time tomography system, the tomography system comprising:

[0018] An anisotropic diffusion equation construction module is used to process the ground seismic profile to obtain ground seismic structure information data, and construct an anisotropic diffusion equation based on the ground seismic structure information data;

[0019] An initial VSP velocity model solving module is used to solve the anisotropic diffusion equation using one-dimensional velocity data as input to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data;

[0020] A time residual matrix calculation module is used to extract the first arrival travel time measurement value in the VSP seismic data; calculate the first arrival travel time calculation value by using the fast scanning method according to the initial VSP velocity model; subtract the first arrival travel time measurement value from the first arrival travel time calculation value to obtain a time residual matrix;

[0021] The iterative update module is used to solve the inversion equation with the time residual matrix as input to obtain an updated VSP velocity model; determine whether the updated VSP velocity model meets the predetermined accuracy; if so, use the updated VSP velocity model as the final VSP velocity model; if not, use the updated VSP velocity model as the initial VSP velocity model in the next iteration, and return to the step of "calculating the first arrival travel time calculated value by the fast scanning method according to the initial VSP velocity model".

[0022] According to the specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0023] This embodiment is used to provide a VSP first arrival traveltime tomography method and system, firstly, the ground seismic profile is processed to obtain ground seismic structure information data, and an anisotropic diffusion equation is constructed according to the ground seismic structure information data, and then the one-dimensional velocity data is used as input to solve the anisotropic diffusion equation to obtain an initial VSP velocity model. Then, the first arrival traveltime measurement value in the VSP seismic data is extracted, and the first arrival traveltime calculation value is calculated by using a fast scanning method according to the initial VSP velocity model, and the first arrival traveltime measurement value and the first arrival traveltime calculation value are subtracted to obtain a time residual matrix, and the time residual matrix is ​​used as input to solve the inversion equation to obtain an updated VSP velocity model, until the updated VSP velocity model meets the predetermined accuracy, and the updated VSP velocity model is used as the final VSP velocity model, so that the ground seismic structure information data is used as a priori knowledge to constrain the inversion process of the VSP velocity model, and improve the accuracy of the VSP first arrival traveltime tomography. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0025] Figure 1 A method flow chart of the tomography method provided in Example 1 of the present invention;

[0026] Figure 2A schematic diagram of the tomography method provided in Example 1 of the present invention;

[0027] Figure 3 A schematic diagram of a two-dimensional layered model provided in Example 1 of the present invention;

[0028] Figure 4 A schematic diagram of superimposed records of a two-dimensional layered model provided in Example 1 of the present invention;

[0029] Figure 5 A schematic diagram of VSP seismic records and initial arrival trip picking provided in Example 1 of the present invention;

[0030] Figure 6 This is a schematic diagram of logging velocity for wellbore position extraction provided by Example 1 of the present invention, wherein: Figure 6 (a) is the true velocity from the model, Figure 6 (b) is the acoustic logging velocity obtained after adding random noise, Figure 6 (c) From Figure 6 The smoothed velocity obtained in (b);

[0031] Figure 7 A schematic diagram of an initial model constructed by diffusing the smoothed logging velocity along the structural direction provided in Example 1 of the present invention;

[0032] Figure 8 A schematic diagram of a velocity model obtained by tomographic inversion provided in Example 1 of the present invention;

[0033] Fig. 9 A schematic diagram of velocity update for tomographic inversion provided in Example 1 of the present invention;

[0034] Fig.10 This is a system block diagram of a tomography system provided in Embodiment 2 of the present invention. DETAILED DESCRIPTION

[0035] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only 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 creative work are within the scope of protection of the present invention.

[0036] The purpose of the present invention is to provide a VSP first-arrival traveltime tomography method and system, which utilizes ground seismic structure information data as prior knowledge to constrain the inversion process of the VSP velocity model and improve the accuracy of VSP first-arrival traveltime tomography.

[0037] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] Embodiment 1:

[0039] After years of development, the application research of VSP technology has been continuously deepened, from zero offset VSP technology to variable offset VSP technology, three-dimensional VSP technology, multi-component VSP technology, inverse VSP technology and Walkaway VSP technology, and with the development of DAS technology, VSP technology has been more widely promoted and applied. However, the amount of VSP seismic data is less than that of ground seismic data, the range of detection points is small, and the target imaging range is limited to the well area. Therefore, the comprehensive use of prior information constraints to improve the inversion accuracy is still the key direction of VSP seismic data inversion research, and it is also a difficult problem that needs to be solved urgently in production practice.

[0040] VSP seismic data has the advantages of high precision and high resolution, but the exploration range of VSP is limited, and it is difficult to accurately invert the medium velocity model. Although the accuracy of ground seismic data is not as good as that of VSP seismic data, it has a wide exploration range. Therefore, the use of ground seismic data and VSP seismic data for joint tomography can achieve complementary advantages among various data, thereby obtaining high-quality inversion results. At present, the research on the joint inversion method of surface seismic data and VSP seismic data includes: (1) using the tomographic information of well logging, VSP and surface seismic data travel time to invert the seismic velocity model; (2) offsetting the surface seismic data, VSP seismic data and cross-hole seismic records to the same imaging space to perform stereo seismic imaging and obtain high-quality imaging of the area near the wellbore; (3) using VSP seismic data and surface seismic data to perform joint tomographic inversion to obtain the velocity of the underground medium; (4) using the first arrival travel time of VSP seismic data to invert anisotropic parameters to assist in establishing a surface seismic anisotropic medium model; (5) carrying out joint imaging of cross-hole seismic, VSP and surface seismic data, performing pre-stack depth migration on seismic data of different scales, and applying the adaptive stacking method after normalization to achieve joint imaging between data of different scales. Many joint inversion methods attempt to solve the shortcomings of VSP seismic data and surface seismic data, and improve the accuracy of VSP seismic data modeling by combining surface seismic data. This embodiment is still based on the idea of ​​comprehensively utilizing ground seismic and VSP seismic data, and proposes a new VSP first-arrival traveltime tomography method under the constraints of ground seismic structural information data. It mainly conducts in-depth research on the use of anisotropic diffusion equation modeling and tomography precondition regularization, so as to achieve the purpose of using prior ground seismic structural information data to improve the accuracy of VSP tomography.

[0041] This embodiment is used to provide a VSP first arrival travel time tomography method under the constraint of ground seismic structure information data, and improves the accuracy of VSP seismic data tomography by introducing prior ground seismic structure information data, such as Figure 1 and Figure 2 As shown, the tomography method comprises:

[0042] S1: Processing the ground seismic profile to obtain ground seismic structure information data, and constructing an anisotropic diffusion equation according to the ground seismic structure information data;

[0043] The ground seismic profile of this embodiment includes a time domain ground seismic profile and a depth domain ground seismic profile. The ground seismic profile is processed to obtain ground seismic structure information data, which may include:

[0044] (1) Determine whether the ground seismic profile is a time domain ground seismic profile;

[0045] (2) If yes, first convert the time domain ground seismic profile into a depth domain ground seismic profile, then perform dip scanning on the depth domain ground seismic profile to obtain ground seismic structural information data;

[0046] (3) If not, that is, the surface seismic profile is a depth-domain surface seismic profile, then the depth-domain surface seismic profile is directly subjected to an inclination scan to obtain surface seismic structural information data.

[0047] This embodiment performs an inclination scan on the depth domain ground seismic profile to obtain ground seismic structural information data. If the ground seismic profile is a time domain ground seismic profile, it needs to be converted into a depth domain ground seismic profile first. The ground seismic structural information data refers to the structural information of the ground seismic profile, which is used to constrain the tomography process. The ground seismic structural information data uses the standard SEGY format. The file in the standard SEGY format generally includes three parts. The first part is the EBCDIC file header, which is used to save some information describing the seismic data body; the second part is the binary file header, which is used to store key information describing the SEGY file, such as data format, number of sampling points, sampling interval and other information; the third part is the actual seismic trace, and each seismic trace contains trace header information and seismic trace data. Alternatively, the ground seismic structural information data uses other data formats containing the position coordinates of the longitudinal survey line (Xline) and the transverse survey line (Inline), the longitudinal sampling interval and the sampling point information.

[0048] The relationship between the VSP seismic data and the ground seismic data used in this embodiment is that the ground seismic imaging profile survey line contains or is adjacent to the VSP well position. The structural information of the ground seismic profile can be used to improve the inversion accuracy during the tomographic imaging of the VSP seismic data travel time. Therefore, the smoothing step cannot destroy the structural information, that is, a smoothing method that protects the ground seismic structural information must be used. The anisotropic diffusion equation based on partial differential equations is a type of adaptive smoothing method. Different smoothing methods can be adopted according to different image contents, and the solution of the equation is the smoothing result.

[0049] The earliest anisotropic diffusion equation can be traced back to the PM model proposed by Perona and Malik, which is in the following form:

[0050]

[0051] Wherein, u(x, z; t) is the seismic image at the position point (x, z) at time t; div is the divergence operator; g(·) is the non-negative decreasing diffusion coefficient; is the gradient operator; u0(x,z) is the original seismic image.

[0052] Weicket analyzed the structure of the image data, converted the diffusion coefficient into a diffusion tensor related to directional characteristics, and constructed the following tensor diffusion equation model:

[0053]

[0054] The above formula is the anisotropic diffusion equation used in this embodiment. Among them, u(x, z; t) is the seismic image at the position point (x, z) at time t obtained from the one-dimensional velocity data; div is the divergence operator; D is the diffusion tensor, which can also be called the tensor diffusion coefficient, and its elements are based on the structure tensor S ρ The local structure information of the extracted image is used to control the diffusion direction and the diffusion rate in the corresponding direction; S ρ is the structure tensor; is the gradient operator; u0(x, z) is the initial seismic image at the (x, z) position point obtained from the one-dimensional velocity data.

[0055] This embodiment uses the above-obtained ground seismic structure information data to solve the diffusion tensor to construct a complete anisotropic diffusion equation. Then, constructing the anisotropic diffusion equation according to the ground seismic structure information data may include:

[0056] (1) Calculate the structural tensor of ground seismic tectonic information data;

[0057] u σ =u*Gσ , where S ρ G is the structural tensor of ground seismic structural information data; ρ represents the Gaussion function with scale ρ; is the gradient operator; u is the seismic image at time t at the position (x, y) obtained from the ground seismic structure information data; G σ represents a Gaussion function with a scale of σ, which is used to avoid the influence of noise when calculating the gradient; G ρ With the squared gradient tensor The purpose of convolution is to take surrounding information into account and avoid canceling out each other when determining edge directions due to the same direction and opposite signs.

[0058] The structure tensor makes each point in the two-dimensional image correspond to a semi-positive definite matrix related to the first-order partial differential information, then the structure tensor is:

[0059]

[0060] Among them, S ρ is the structure tensor; G ρ is a Gaussion function with a scale of ρ; is the gradient operator; u σ =u*G σ , u is the seismic image at the position (x, y) obtained from the ground seismic structure information data; G σ is a Gaussion function with a scale of σ.

[0061] (2) According to the decomposition theorem of matrix eigenvalues, the structure tensor is decomposed to obtain the first eigenvalue and the first eigenvector of the structure tensor;

[0062] After eigenvalue decomposition, the structure tensor can be expressed as:

[0063]

[0064] Among them, v1, v2 are the first eigenvectors; λ1, λ2 are the first eigenvalues ​​corresponding to the first eigenvectors; v1 represents parallel The direction of the maximum rate of change is v2 means perpendicular to The direction of the smallest rate of change is

[0065] The specific solution form of the first eigenvector is as follows:

[0066]

[0067] The corresponding eigenvalues ​​are:

[0068]

[0069] (3) The first eigenvalue is normalized and regularized to obtain a second eigenvalue; the second eigenvalue and the first eigenvector are used as the eigenvalue and eigenvector of the diffusion tensor, the diffusion tensor in the anisotropic diffusion equation is constructed, and the anisotropic diffusion equation is obtained.

[0070] The eigenvector of the diffusion tensor is the same as that of the structure tensor, but the eigenvalue is different. The eigenvalue needs to be calculated according to the specific requirements of image processing. Weickert gave the eigenvalue of the diffusion tensor based on the requirements of edge enhancement and coherence enhancement. The eigenvalues ​​u1 and u2 independently control the diffusion behavior in the v1 and v2 directions respectively. The selected eigenvalues ​​are:

[0071]

[0072] The diffusion tensor obtained by the above formula uses the eigenvector of the structure tensor to accurately orient the changes in the image. At the same time, the size of the change along the feature direction (i.e., the size of the eigenvalue) can be defined according to the needs of image processing, and no longer depends solely on the local structure tensor or local structure of the image.

[0073] Different from image processing, seismic exploration problems need to meet the following two requirements: the diffusion tensor is infinitely smooth when viewed as a function of the structure tensor, and the eigenvalue of the diffusion tensor has a positive lower bound. Based on the two requirements and combined with the actual problems of seismic exploration, this embodiment solves the eigenvalue of the diffusion tensor through the following two steps:

[0074] (3.1) Perform L2 norm normalization on the first eigenvalue to obtain the normalized eigenvalue;

[0075]

[0076] Among them, A and B are normalized eigenvalues, and λ=(λ1, λ2) is the eigenvalue vector.

[0077] (3.2) The normalized eigenvalue is regularized according to the lower bound of the diffusion tensor eigenvalue to obtain the second eigenvalue.

[0078] Considering that the diffusion tensor eigenvalue has a positive lower bound, the normalized eigenvalue is normalized to the range [α, 1] by the following formula:

[0079]

[0080] According to the relationship between the diffusion tensor and the structure tensor, the regular diffusion tensor can be constructed using the new second eigenvalue and the eigenvector of the structure tensor. After obtaining the regular diffusion tensor, the anisotropic diffusion equation can be established.

[0081] The small observation range of VSP seismic exploration leads to the incompleteness of VSP data and the ill-posedness of tomography. Therefore, it is very necessary to add prior information to tomography. By constructing the diffusion tensor of the anisotropic diffusion equation, the ground seismic structure information is introduced into the VSP seismic data tomography, and the tomography is constrained by the prior information.

[0082] S2: using one-dimensional velocity data as input, solving the anisotropic diffusion equation to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data;

[0083] The VSP velocity data and the logging velocity data are uniformly referred to as one-dimensional velocity data, and the one-dimensional velocity data is in text format, including the well name, location and depth, and velocity information. The one-dimensional velocity data (VSP velocity data or logging velocity data) is checked, and unreasonable outliers (negative values, etc.) are removed to obtain reasonable one-dimensional velocity data, which is the one-dimensional velocity data in S2 of this embodiment. The specific solution method for solving the anisotropic diffusion equation with the one-dimensional velocity data as input is as follows: construct an empty matrix, fill the one-dimensional velocity data into the corresponding position in the matrix according to the position information of the one-dimensional velocity data, and use the filled matrix as the initial image u0 in the anisotropic diffusion equation, and then solve the equation for diffusion to obtain the initial VSP velocity model.

[0084] In this embodiment, ground seismic structure information data is obtained through the ground seismic profile in the depth domain, and the ground seismic structure information data is used to construct the diffusion tensor in the anisotropic diffusion equation to construct the anisotropic diffusion equation, and then the anisotropic diffusion equation is solved based on the one-dimensional velocity data, and the velocity at the wellbore is diffused into a two-dimensional or three-dimensional initial VSP velocity model. After the initial VSP velocity model is obtained, it is used as input for VSP first arrival travel time tomography. Based on the one-dimensional velocity data, the initial VSP velocity model is established by diffusion along the ground seismic large-scale structure direction using the anisotropic diffusion equation, and the ground seismic structure information data can be used to provide a good initial model for tomography.

[0085] S3: Extract first arrival travel time measurements from VSP seismic data;

[0086] The VSP seismic data of this embodiment uses a standard SEGY format, or the VSP seismic data uses a data format that includes the travel time of the first arrival wave of each data, the shot number, the coordinates of the source point and the detection point, the time sampling interval and the number of sampling points. The positions of the source point and the detection point can be used to determine the positions of seismic wave excitation and reception. The travel time of the first arrival wave of the VSP seismic data is subtracted from the travel time obtained by the fast scanning method to obtain the travel time error, and the inversion equation is further solved.

[0087] S4: according to the initial VSP velocity model, using a fast scanning method to calculate the first arrival travel time;

[0088] Specifically, S4 may include: constructing an eikonal equation according to the initial VSP velocity model; solving the eikonal equation using a fast scanning method to obtain a calculated value of the first arrival travel time.

[0089] More specifically, under the assumption of "high frequency approximation" (i.e., the wavelength of seismic waves is much smaller than the inhomogeneity of the medium, and the seismic wave energy is mainly propagated along rays), the elastic wave equation can be simplified to the Eikonal equation, which gives the relationship between the magnitude of the travel time gradient and the slowness (the inverse of the velocity) at any point. The expression of the Eikonal equation in two-dimensional space is:

[0090]

[0091] in, is the gradient operator, T is the travel time, T(x) is the travel time of the point x, and c(x) is the propagation velocity of the seismic wave at the point x in the medium. The travel time gradient and slowness of the direct wave satisfy the expression of the eikonal equation, and at the source point, T(x s )=0,x s is the location of the earthquake source,

[0092] The fast scanning method divides the propagation direction of the travel time field into a finite number of groups based on the causal relationship. For each group, the Gauss-Seidel iterative equation is used to solve the discretized discrete equation of the upwind difference format. The discrete equation of the upwind difference format is:

[0093]

[0094] Among them, t i,j is the travel time at the grid point (i, j); i, j are the number of grid points in the x and z directions respectively, s i,j is the slowness (inverse of the velocity) at the grid point (i, j), h is the spatial step size, and

[0095]

[0096] Among them, f refers to () + If the data or expression in f is greater than 0, then (f) + =f, if f is less than or equal to 0, then (f) + =0.

[0097] In two-dimensional space, the scanning directions can be divided into four groups (upper right, upper left, lower left and lower right). The first arrival information is propagated along the four directions. The propagation travel time of the seismic wave is solved by four groups of Gauss-Seidel iterations in different orders combined with the upwind difference format.

[0098] The steps for calculating the travel time of the earthquake first arrival wave by the fast scanning method are:

[0099] (1) Initialization: Assign an initial value of 0 to the focal point and assign precise values ​​to the four points around the focal point to avoid first-order numerical errors. The four points around the focal point refer to the four points above, below, left and right of the focal point. The precise value of the travel time can be obtained by solving the Eikonal equation. The initial values ​​of the grid points other than the focal point are set to a large value, which is much larger than the travel time that can be calculated at all calculation points. The travel time at these points will be updated in subsequent calculations.

[0100] (2) Iteration: In four directions, the Gauss-Seidel type iteration method is used to solve the discrete nonlinear equations. At each grid point where the travel time is not fixed, the current travel time is calculated and compared with the original value, and the smaller value is taken as the current value;

[0101] (3) Convergence judgment: If the convergence condition || t is met new -t old ||≤δ, the calculation ends; if it is not satisfied, the iteration continues until convergence. The current travel time of each grid point other than the source point and the above four surrounding points is calculated using the discrete equation of the upwind difference format.

[0102] S5: subtracting the first arrival travel time measurement value from the first arrival travel time calculation value to obtain a time residual matrix; using the time residual matrix as input, solving the inversion equation to obtain an updated VSP velocity model;

[0103] The basic principle of tomography is to first give an initial slowness model S (the inverse of the velocity model), construct a ray path matrix L based on ray path tracing, solve the eikonal equation by the fast scanning method to obtain the calculated value of the first arrival travel time of the seismic wave, and subtract the first arrival travel time measurement value from the first arrival travel time calculated value to obtain the travel time residual matrix ΔT, thereby forming the inversion equation:

[0104] LΔS=ΔT;

[0105] Among them, L is the ray path matrix; ΔS is the slowness model, slowness is the inverse of speed; ΔT is the time residual matrix.

[0106] Solving the above large sparse equations can obtain the medium slowness, that is, the updated VSP velocity model. This equation is an ill-conditioned large coefficient linear algebraic equation system, which can be solved by algebraic reconstruction technology (ART), instantaneous iterative reconstruction technology (SIRT), damped least squares method (LSQR) and other methods.

[0107] S6: Determine whether the updated VSP velocity model meets the predetermined accuracy; if so, use the updated VSP velocity model as the final VSP velocity model; if not, use the updated VSP velocity model as the initial VSP velocity model in the next iteration, and return to the step of "calculating the first arrival travel time by the fast scanning method according to the initial VSP velocity model" to obtain the final inversion velocity through continuous iteration.

[0108] As an optional implementation, after obtaining the updated VSP velocity model and before determining whether the updated VSP velocity model satisfies a predetermined accuracy, the tomography method of this embodiment further includes: determining whether the number of iterations reaches a predetermined number; if so, performing preconditional smoothing on the slowness model in the inversion equation using the anisotropic diffusion equation to obtain a new inversion equation, and using the new inversion equation as the inversion equation in the next iteration, while setting the number of iterations to 0; if not, not changing the inversion equation and continuing the iteration.

[0109] Specifically, during the model updating process, the anisotropic diffusion equation can be constructed using the ground seismic structure information data to perform precondition smoothing on the slowness model:

[0110] ΔS′=QΔS;

[0111] The above formula is a schematic representation of the anisotropic diffusion equation, which is the filter that realizes the above process. Among them, Q is the preconditioning operator, and QΔS means that the model ΔS is smoothly diffused using the anisotropic diffusion equation. In order to improve the inversion efficiency, it is not necessary to smooth the model through the anisotropic diffusion equation every time it is updated. The model can be intermittently smoothed along the ground seismic structure direction before entering the next model update.

[0112] This embodiment uses pre-conditioned regularization to smoothly diffuse the model along the seismic structure direction, which can introduce ground seismic structure information to improve the quality of tomography. At the same time, the velocity model is intermittently diffused and smoothed along the seismic structure direction during the iteration process to more efficiently construct a high-precision VSP velocity model.

[0113] In this embodiment, a dip scan is first performed on a surface seismic profile in the depth domain to obtain surface seismic structure information. Based on one-dimensional VSP velocity data or logging velocity data, an initial velocity model is constructed along the surface seismic large-scale structure direction using anisotropic diffusion equations. The rapid scanning method is used to solve the eikonal equation to obtain the earthquake first arrival travel time. A large sparse matrix linear algebraic equation group is solved to invert the velocity model, thereby achieving efficient and accurate VSP data tomography. The surface seismic structure information is used to perform pre-conditional tomography. In the iterative step of tomography inversion, the velocity model is intermittently diffused and smoothed along the surface seismic structure direction to improve the accuracy of VSP data tomography, thereby using seismic prior structure information to improve the accuracy of VSP first arrival travel time tomography.

[0114] Next, we use a two-dimensional layered model to test the validity of the first arrival traveltime tomography method for preconditioned VSP data. Figure 3 As shown in the figure, the size is 9000m×6000m, the density is constant, and the velocity contains multiple layers of interface constraints with undulating changes. The model velocities are 2500m / s, 3000m / s, 3400m / s, 4100m / s and 4500m / s from top to bottom. The model structure is complex, the stratum undulation varies greatly, and it is also difficult to perform tomography. VSP seismic data is obtained through the velocity model ( Figure 3 ) obtained by forward modeling, a total of 121 shots, with a shot spacing of 50 meters, a source burial depth of 0.5 meters, and evenly distributed within the range of 0.6 kilometers to 6.6 kilometers. The wellbore is located 0.5 kilometers in the horizontal direction, and the receiving points are evenly distributed within the depth range of 1 kilometer to 5 kilometers.

[0115] Figure 4 It is a simulated stacked profile, which is mainly used to generate underground structural information in the area where the velocity model is located. The surface seismic profile is a depth domain seismic profile. In actual production, most seismic profiles are pre-stack time migration profiles. The corresponding depth domain profile can be obtained through time-to-depth conversion. Figure 5 It is a VSP seismic record obtained based on the velocity model. In order to perform first-arrival traveltime tomography of VSP seismic data, it is necessary to pick up the first-arrival traveltime. The circle in the figure represents the picked first-arrival arrival time, which serves as the input data for first-arrival traveltime tomography.

[0116] In order to simulate the steps of real processing, Figure 6 Indicates the location of the VSP well, Figure 6 (a) is the true velocity from the model, Figure 6 (b) is the acoustic logging velocity obtained after adding random noise, Figure 6 (c) is from Figure 6 (b) The smoothed velocity obtained in this embodiment is based on the one-dimensional velocity model and uses the anisotropic diffusion equation to obtain the initial model of tomography, such as Figure 7 As shown in the figure, the initial velocity model is formed by using the structure of the seismic stacked section as a constraint diffusion, which obviously contains ground structure information and is generally closer to the real velocity field. Figure 8 It is the velocity model obtained by tomographic inversion based on the initial model. Fig. 9 is the update amount of tomographic inversion. The tomographic update results are significantly different from the initial model, mainly because the update amount is larger at the structural boundary and there are structural traces of superimposed records in the update amount.

[0117] Embodiment 2:

[0118] This embodiment is used to provide a VSP first arrival travel time tomography system, such as Fig.10 As shown, the tomography system comprises:

[0119] Anisotropic diffusion equation construction module M1, used for processing ground seismic profiles to obtain ground seismic structure information data, and constructing anisotropic diffusion equations according to the ground seismic structure information data;

[0120] The initial VSP velocity model solving module M2 is used to solve the anisotropic diffusion equation using one-dimensional velocity data as input to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data;

[0121] The time residual matrix calculation module M3 is used to extract the first arrival travel time measurement value in the VSP seismic data; calculate the first arrival travel time calculation value by using the fast scanning method according to the initial VSP velocity model; subtract the first arrival travel time measurement value from the first arrival travel time calculation value to obtain the time residual matrix;

[0122] The iterative update module M4 is used to solve the inversion equation with the time residual matrix as input to obtain an updated VSP velocity model; determine whether the updated VSP velocity model meets the predetermined accuracy; if so, use the updated VSP velocity model as the final VSP velocity model; if not, use the updated VSP velocity model as the initial VSP velocity model in the next iteration, and return to the step of "calculating the first arrival travel time calculated value by the fast scanning method according to the initial VSP velocity model".

[0123] Each embodiment in this specification focuses on the differences from other embodiments, and the same or similar parts between the embodiments can be referred to each other. For the system disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0124] The principles and implementation methods of the present invention are described in this article using specific examples. The description of the above embodiments is only used to help understand the method and core idea of ​​the present invention. At the same time, for those skilled in the art, according to the idea of ​​the present invention, there will be changes in the specific implementation methods and application scope. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A VSP first arrival travel time tomography method, characterized in that: The tomography method comprises: Processing the ground seismic profile to obtain ground seismic structure information data, and constructing an anisotropic diffusion equation based on the ground seismic structure information data; The anisotropic diffusion equation is solved using one-dimensional velocity data as input to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data; Extract first arrival traveltime measurements from VSP seismic data; According to the initial VSP velocity model, the first arrival travel time calculation value is calculated by using the fast scanning method; Subtracting the first arrival travel time measurement value from the first arrival travel time calculation value to obtain a time residual matrix; Taking the time residual matrix as input, solving the inversion equation to obtain an updated VSP velocity model; Determining whether the updated VSP velocity model meets a predetermined accuracy; If yes, the updated VSP velocity model is used as the final VSP velocity model; If not, the updated VSP velocity model is used as the initial VSP velocity model in the next iteration, and the process returns to the step of "calculating the first arrival travel time by using the fast scanning method according to the initial VSP velocity model"; The constructing of an anisotropic diffusion equation according to the ground seismic structure information data specifically includes: Calculating the structural tensor of the ground seismic structure information data; Performing eigenvalue decomposition on the structure tensor to obtain a first eigenvalue and a first eigenvector of the structure tensor; Normalizing and regularizing the first eigenvalue to obtain a second eigenvalue; The second eigenvalue and the first eigenvector are used as the eigenvalue and eigenvector of the diffusion tensor to construct the diffusion tensor in the anisotropic diffusion equation, thereby obtaining the anisotropic diffusion equation.

2. The tomography method according to claim 1, characterized in that: The processing of the ground seismic profile to obtain ground seismic structure information data specifically includes: Determining whether the ground seismic profile is a time-domain ground seismic profile; If yes, converting the time-domain ground seismic profile into a depth-domain ground seismic profile; The depth domain ground seismic profile is subjected to an inclination scan to obtain ground seismic structure information data.

3. The tomographic imaging method according to claim 1, characterized in that: The structure tensor is: Among them, S ρ is the structure tensor; G ρ is a Gaussian function with scale ρ; is the gradient operator; u σ =u*G σ , u is the seismic image at the position (x, y) obtained from the ground seismic structure information data; G σ is a Gaussian function with scale σ.

4. The tomography method according to claim 1, characterized in that: The normalizing and regularizing the first eigenvalue to obtain the second eigenvalue specifically includes: Performing L2 norm normalization on the first eigenvalue to obtain a normalized eigenvalue; The normalized eigenvalue is regularized according to a lower bound of the diffusion tensor eigenvalue to obtain a second eigenvalue.

5. The tomography method according to claim 1, characterized in that: The anisotropic diffusion equation is: Where u(x, z; t) is the seismic image at the position point (x, z) at time t obtained from the one-dimensional velocity data; div is the divergence operator; D is the diffusion tensor; S ρ is the structure tensor; is the gradient operator; u0(x, z) is the initial seismic image at the (x, z) position point obtained from the one-dimensional velocity data.

6. The tomography method according to claim 1, characterized in that: The method of calculating the initial arrival travel time using the rapid scanning method according to the initial VSP velocity model specifically includes: constructing an eicoplanin equation according to the initial VSP velocity model; The rapid scanning method is used to solve the eikonal equation to obtain the calculated value of the first arrival travel time.

7. The tomography method according to claim 1, characterized in that: The inversion equation is: LΔS=ΔT Among them, L is the ray path matrix; ΔS is the slowness model, slowness is the inverse of speed; ΔT is the time residual matrix.

8. The tomography method according to claim 1, characterized in that: After obtaining the updated VSP velocity model and before determining whether the updated VSP velocity model satisfies a predetermined accuracy, the tomography method further includes: Determine whether the number of iterations reaches a predetermined number; If so, the anisotropic diffusion equation is used to perform precondition smoothing on the slowness model in the inversion equation to obtain a new inversion equation, and the new inversion equation is used as the inversion equation in the next iteration; and the number of iterations is set to 0.

9. A VSP first arrival travel time tomography system, characterized in that: The tomography system comprises: The anisotropic diffusion equation construction module is used to process the ground seismic profile to obtain ground seismic structure information data, and construct the anisotropic diffusion equation according to the ground seismic structure information data; the construction of the anisotropic diffusion equation according to the ground seismic structure information data specifically includes: Calculating the structural tensor of the ground seismic structure information data; Performing eigenvalue decomposition on the structure tensor to obtain a first eigenvalue and a first eigenvector of the structure tensor; Normalizing and regularizing the first eigenvalue to obtain a second eigenvalue; Using the second eigenvalue and the first eigenvector as the eigenvalue and eigenvector of the diffusion tensor, constructing the diffusion tensor in the anisotropic diffusion equation, and obtaining the anisotropic diffusion equation; An initial VSP velocity model solving module is used to solve the anisotropic diffusion equation using one-dimensional velocity data as input to obtain an initial VSP velocity model; the one-dimensional velocity data is VSP velocity data or logging velocity data; A time residual matrix calculation module is used to extract the first arrival travel time measurement value in the VSP seismic data; calculate the first arrival travel time calculation value by using the fast scanning method according to the initial VSP velocity model; subtract the first arrival travel time measurement value from the first arrival travel time calculation value to obtain a time residual matrix; The iterative update module is used to solve the inversion equation with the time residual matrix as input to obtain an updated VSP velocity model; determine whether the updated VSP velocity model meets the predetermined accuracy; if so, use the updated VSP velocity model as the final VSP velocity model; if not, use the updated VSP velocity model as the initial VSP velocity model in the next iteration, and return to the step of "calculating the first arrival travel time calculated value by using the fast scanning method according to the initial VSP velocity model".

Images