Velocity model smoothing method based on least squares convolution fitting

Abstract

The present invention relates to a velocity model smoothing processing method based on least squares convolution fitting. It mainly solves the problem that the conventional commercial software velocity model smoothing method for unconventional oil and gas exploration targets will cause velocity model distortion, resulting in imaging and structural accuracy that cannot meet exploration needs. The velocity model smoothing processing method disclosed in the present invention determines the smoothing parameters S in three directions through the target layer depth H (or time domain T and target layer average velocity V) X , S Y , S Z , and then the least square convolution fitting algorithm is used to smooth the three-dimensional velocity model. The posterior effect is good. The smoothed velocity model is stable and consistent with the structural trend. The migration imaging is accurate, and the gather event axes are flattened. Using the method of the present invention, the input velocity model observation information can be better retained, and the final imaging and structural accuracy are effectively improved, providing data support for unconventional oil and gas exploration.

Description

Technical field:

[0001] The present invention relates to the technical field of seismic data processing, and in particular to a velocity model smoothing processing method based on least square convolution fitting. The method provides an accurate and smooth velocity model for PSDM or PSTM. Background technology:

[0002] Before the 21st century, exploration targets were mainly structural exploration, and the imaging accuracy requirements were relatively low. Therefore, conventional prestack time migration did not have high requirements for velocity fields. Since the 21st century, exploration targets have gradually shifted from structural exploration to unconventional oil and gas exploration, which puts forward higher requirements for imaging accuracy. The most important factor affecting the accuracy of seismic migration imaging is velocity. Whether it is prestack time migration PSTM or prestack depth migration PSDM, the accuracy and rationality of the migration velocity directly determine the accuracy of the migration results. However, the interactive velocity analysis of mainstream commercial software for seismic data processing currently completes the analysis at a certain distance interval (integer multiples of the trace spacing). The obtained time / depth-velocity pairs need to be interpolated and smoothed before they can be used for overall migration processing. When the velocity smoothing result changes dramatically laterally, the integral migration method will be unstable. In mainstream commercial seismic data processing software, velocity smoothing algorithms and implementation methods are not taken as important research and development objects. Commonly used smoothing algorithms such as median filtering, Laplace filtering, and sliding smoothing algorithms are all commonly used smoothing methods in mainstream seismic data processing software. The default smoothing parameters usually do not consider the complexity of the actual processing object structure. Due to the limitation of the smoothing method, it is easy to cause abnormal local changes in the smoothed velocity model or damage the effective information. For seismic processing, the smoothing of the migration velocity field has not attracted enough attention from practitioners, but in actual production, the rationality of the velocity field smoothing result directly determines the success or failure of the whole processing flow. Summary of the invention:

[0003] The present invention aims at the current unconventional oil and gas exploration targets. Conventional commercial software velocity model smoothing methods may cause local velocity anomalies caused by insufficient velocity model smoothing or velocity model distortion caused by excessive velocity model smoothing, thereby resulting in the problem that imaging and structural accuracy cannot meet exploration needs. A velocity model smoothing processing method based on least squares convolution fitting is provided. Based on the technology of the present invention, the smoothed velocity model is stable, the input velocity model observation information can be better retained, the velocity model smoothing accuracy is improved, and finally the imaging and structural accuracy are effectively improved, providing data support for unconventional oil and gas exploration.

[0004] The present invention solves the problem by the following technical solution: the velocity model smoothing processing method based on least squares convolution fitting comprises the following steps:

[0005] S1. Select the target block and select the migration method based on the geological requirements of the target block;

[0006] S2. Based on the selected migration method, the optimal sliding window parameter S relative to the target layer is calculated through the initial 3D velocity model. X , S Y , S Z ; The optimal sliding window parameter S is obtained X , S Y , S Z , determine S X ×S Y ×S Z Sliding window;

[0007] Where: S X is the smoothing window parameter of the main survey line direction of the 3D velocity model;

[0008] S Y is the smoothing window parameter of the contact line direction of the 3D velocity model;

[0009] S Z is the vertical smoothing window parameter of the three-dimensional velocity model;

[0010] S3. The S calculated in step S2 is X , S Y , S Z , for S in the sliding window X ×S Y ×S Z The speed values ​​of the sample points are calculated by multivariate n-order polynomial fitting, and the S corresponding to the polynomial is determined. X ×S Y ×S Z Weight coefficient matrix A;

[0011] S4 takes the S determined in step S3 X ×S Y ×S Z The weight matrix A is applied to S3 X ×S Y ×S Z All elements of the sliding window can obtain the best fitting value at the center point D(i,j,k) of the sliding window;

[0012] S5. Smoothing the entire initial three-dimensional velocity model along the X, Y, and Z directions of the initial three-dimensional velocity model respectively. After the smoothing process is completed, a smoothed velocity model is obtained; and the smoothed velocity model is output for migration processing.

[0013] Preferably, the migration method comprises: using a prestack depth migration processing method to smooth the depth domain layer velocity model and using a prestack time migration processing method to smooth the time domain migration velocity model.

[0014] Preferably, the method of smoothing the depth domain layer velocity model by using the prestack depth migration processing method comprises:

[0015] First, determine the depth H of the imaging target layer;

[0016] Calculate the smoothing window parameter S of the main survey line direction of the 3D velocity model X ;

[0017] And the smoothing window number S of the contact line direction of the 3D velocity model Y ;

[0018] and the vertical smoothing window parameter S of the 3D velocity model Z ;

[0019] Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness;

[0020] and / or,

[0021] S X =4H / Δx+1;

[0022] S Y =4H / Δy+1;

[0023] S X , S Y , S Z All are odd numbers;

[0024] Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in units of m;

[0025] H is the depth of the imaging target layer, in meters.

[0026] Preferably, the vertical smoothing parameter S of the three-dimensional velocity model is Z , determined according to the requirements for vertical smoothness of the model, the vertical smoothing parameter S Z Take S X and S Y The larger value.

[0027] Preferably, the method of smoothing the time domain migration velocity model by using the prestack time migration processing method comprises:

[0028] First, determine the imaging target layer time T and the average velocity V of the target layer;

[0029] Calculate the smoothing parameter S of the main survey line direction of the 3D velocity model X ;

[0030] And the smoothing window number S of the contact line direction of the 3D velocity model Y ,

[0031] and the vertical smoothing parameter S of the 3D velocity model Z ;

[0032] Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness;

[0033] and / or,

[0034] S X =2VT / Δx+1;

[0035] S Y =2VT / Δy+1;

[0036] S X , S Y , S Z All are odd numbers;

[0037] Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in units of m;

[0038] T is the imaging target layer time, in seconds;

[0039] V is the average velocity of the target layer, in m / s.

[0040] The S Z The vertical smoothing parameter S is determined according to the requirement for vertical smoothness of the model. Z Take S X and S Y The larger value.

[0041] Preferably, in step S3, X ×S Y ×S Z The method for obtaining the weight coefficient matrix A includes the following steps:

[0042] S71. Determine the least squares fitting polynomial, with S X , S Y , S ZThe number of elements in the x, y, and z directions in the sliding window of size are l, m, and n respectively, and the element at the center is D(i, j, k), and the fitting polynomial can be obtained:

[0043]

[0044] Where:

[0045]

[0046]

[0047]

[0048] P(i',j',k',) is the fitting value;

[0049] (x i' ,y j' ,z k' ) is the position of the element D(i',j',k') in the smoothing window;

[0050] a svt To fit the polynomial coefficients, write it in matrix form: P = AB, where P is the matrix composed of each element in the sliding window after fitting, A is the matrix composed of the least squares fitting polynomial parameters, that is, the weight coefficient matrix, and matrix B is the elements in the sliding window according to the coefficient a in the fitting polynomial. svt The corresponding x i' s y j' t z k' v A matrix of values;

[0051] l, m, and n represent the number of elements in the sliding window in the three directions of x, y, and z respectively; x i' s y j' t z k' v Represents the position of each point in the sliding window (x i' ,y j' ,z k' ) is the value calculated according to the polynomial fitting formula.

[0052] S72. The mean square error E=PD in the sliding window obtained by fitting can be determined; wherein: the matrix E is the mean square error, and the matrix D is a matrix composed of the initial values ​​of each element in the sliding window;

[0053] S73. According to the least squares principle, the root mean square error is minimized, and the weight coefficient matrix A = B can be obtained by transforming the mean square error E. -1D,

[0054] Preferably, in step S3, the polynomial degree n X , S Y , S Z , n is 3.

[0055] Preferably, the method of smoothing the entire three-dimensional velocity data volume along the three directions of X, Y and Z respectively in step S5 comprises the following steps:

[0056] Select the sliding window from the initial position and select the sliding window S in increments of 1 unit. X ×S Y ×S Z Time window, judge sliding S X ×S Y ×S Z Whether the end of the time window reaches the maximum element of the initial three-dimensional velocity model in this direction, if not, repeat the operations of steps S3 and S4; in the X, Y, and Z directions of the initial three-dimensional velocity model, when the sliding S X ×S Y ×S Z When the end of the time window reaches the maximum element of the initial three-dimensional velocity model in this direction, the sliding is ended, and the least square convolution fitting weighted processing of the initial three-dimensional velocity model in this direction is stopped;

[0057] When processing in all three directions stops, the smoothing process of the entire initial three-dimensional velocity data volume is completed.

[0058] The velocity model smoothing processing method based on least squares convolution fitting of the present invention is comparable to the least squares convolution fitting algorithm in terms of noise resistance and anti-outlier capability compared to conventional median filtering, Laplace filtering, and sliding smoothing algorithms, but the least squares convolution fitting method better preserves the observation information of the data. The least squares convolution fitting algorithm is the most commonly used signal processing method and is widely used in fields such as spectral analysis. Least squares convolution fitting (polynomial fitting) is a very mature algorithm, and there are mature programs in software such as MATLAB and various open source databases, which can be quickly implemented through simple variable settings. At the same time, the application object of the present invention is the velocity model in the time domain and the depth domain. The post-stack velocity model is different from the volume of pre-stack seismic data. The post-stack velocity model is very small. Taking a standard segy format velocity model with a recording track length of 10,000 m, 1,000 lines, 1,000 CDP / line, and a vertical sampling interval of 10 m as an example, the model size is less than 10 GB, and there is almost no need to consider the calculation efficiency. The least squares convolution fitting has already realized multi-threaded parallel calculation. The calculation efficiency and the implementation method are not bottleneck problems that restrict the implementation and promotion of the core content of the present invention. Therefore, the least squares convolution fitting algorithm and the implementation process are not elaborated in the present invention. ​

[0059] Compared with the above background technology, the present invention has the following beneficial effects:

[0060] A stable and correct velocity model is the basis of migration imaging. The present invention takes the velocity smoothing algorithm and implementation method as important research and development objects. After smoothing by the method of the present invention, there is no lateral velocity anomaly and vertical velocity reversal in the velocity model profile as a whole, making the velocity field smoothing result more reasonable; making the lateral change trend and vertical change law of the velocity field more consistent with the understanding of geological structure.

[0061] The least squares fitting smoothing method of the present invention is based on the least squares principle and effectively solves the problem encountered in actual production that the local distortion of the velocity model caused by insufficient smoothing or the distortion of the velocity model caused by smoothing overflow leads to the imaging accuracy failing to meet the needs of fine exploration. The smoothed velocity model is stable, the observation information of the input velocity model can be well retained, the degree of consistency with the structural trend is high, the migration imaging accuracy is high, and the structural error is small. In the process of velocity smoothing, the effective information of the velocity pickup of the quality control point can be retained, and the reasonable changes of the structure can be followed to the greatest extent. The smoothed velocity field has a very high degree of consistency with the structure in the horizontal and vertical directions. Under the guidance of the least squares fitting smoothing technology, the fidelity of the amplitude of the final migration result profile is good, the well compliance rate is increased from 0.7 to 0.91, and the imaging accuracy is significantly improved.

[0062] The velocity model smoothed by the method of the present invention has the advantages of simplicity and high efficiency compared with the conventional commercial software velocity model smoothing method. It has been effectively promoted and applied in multiple horizontal well exploration areas in Daqing Oilfield with T2 overlying mudstone interlayer as the exploration target. The velocity model after smoothing by the method of the present invention is stable, the correlation with the structural trend is significantly improved, the accuracy is greatly improved, the final result data imaging has high accuracy and good fidelity, and the well structure error is significantly improved. It has broad application prospects in the field of unconventional oil and gas exploration (tight oil, deep gas, etc.) in Daqing Oilfield.

[0063] The method was promoted and applied in a horizontal well exploration area in Daqing and achieved good results. The final time domain migration results have high fidelity, and the well structure error of depth domain imaging has been significantly improved. Description of the drawings:

[0064] Figure 1 is a flow chart of the velocity model smoothing method of the present invention;

[0065] Figure 2 is the RMS velocity model profile after smoothing by default parameters of commercial seismic data processing software Omega in the embodiment of the present invention;

[0066] Figure 3 is the RMS velocity model profile after smoothing by the method of the present invention in the embodiment of the present invention, Figure 2 Same location;

[0067] Figure 4 In the embodiment of the present invention Figure 3 , Figure 4 The actual seismic profile corresponding to the same position of the velocity profile;

[0068] 5 is a synthetic record of the seismic profile obtained by the conventional velocity smoothing method in an embodiment of the present invention (a), and a synthetic record of the seismic profile obtained by the least squares fitting velocity smoothing method (b). Specific implementation method:

[0069] In order to make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0070] like Figure 1 As shown, a velocity model smoothing processing method based on least squares convolution fitting includes the following steps:

[0071] S1. Select the target block and select the migration method based on the geological requirements of the target block;

[0072] The migration method includes: using a pre-stack depth migration processing method to smooth a depth domain layer velocity model and using a pre-stack time migration processing method to smooth a time domain migration velocity model.

[0073] S2. Based on the selected migration method, the optimal sliding window parameter S relative to the target layer is calculated through the initial 3D velocity model. X , S Y , S Z ; The optimal sliding window parameter S is obtained X , S Y , S Z , determine S X ×S Y ×S Z Sliding window;

[0074] Where: S X is the smoothing window parameter of the main survey line direction of the 3D velocity model;

[0075] S Y is the smoothing window parameter of the contact line direction of the 3D velocity model;

[0076] S Z is the vertical smoothing window parameter of the three-dimensional velocity model;

[0077] S21. The method of smoothing the depth domain layer velocity model by using the prestack depth migration processing method comprises the following steps:

[0078] First, determine the depth H of the imaging target layer;

[0079] Calculate the smoothing window parameter S of the main survey line direction of the 3D velocity model X ;

[0080] And the smoothing window number S of the contact line direction of the 3D velocity model Y ;

[0081] and the vertical smoothing window parameter S of the 3D velocity model Z ;

[0082] The vertical smoothing parameter S of the three-dimensional velocity model Z , determined according to the requirements for vertical smoothness of the model, the vertical smoothing parameter S Z Take S X and S Y The larger value.

[0083] Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness;

[0084] and / or,

[0085] S X =4H / Δx+1;

[0086] S Y =4H / Δy+1;

[0087] S X , S Y , S Z All are odd numbers;

[0088] Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in meters; H is the depth of the imaging target layer, in meters.

[0089] S22. The method of smoothing the time domain migration velocity model using the prestack time migration processing method comprises the following steps:

[0090] First, determine the imaging target layer time T and the average velocity V of the target layer;

[0091] Calculate the smoothing parameter S of the main survey line direction of the 3D velocity model X ;

[0092] And the smoothing window number S of the contact line direction of the 3D velocity modelY ,

[0093] and the vertical smoothing parameter S of the 3D velocity model Z ; S Z The vertical smoothing parameter S is determined according to the requirement for vertical smoothness of the model. Z Take S X and S Y The larger value.

[0094] Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness;

[0095] and / or,

[0096] S X =2VT / Δx+1;

[0097] S Y =2VT / Δy+1;

[0098] S X , S Y , S Z All are odd numbers;

[0099] Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in units of m;

[0100] T is the imaging target layer time, in seconds;

[0101] V is the average velocity of the target layer, in m / s.

[0102] S3. The S calculated in step S2 is X , S Y , S Z , for S in the sliding window X ×S Y ×S Z The velocity values ​​of the sample points are fitted with a multivariate n-order polynomial to determine the S corresponding to the polynomial. X ×S Y ×S Z Weight coefficient matrix A; polynomial degree n X , S Y , S Z , n is 3.

[0103] The S X ×S Y ×S Z ​The method for obtaining the weight coefficient matrix A includes the following steps:

[0104] S31. Determine the least squares fitting polynomial, S X , S Y , S Z The number of elements in the x, y, and z directions in the sliding window of size are l, m, and n respectively, and the element at the center is D(i, j, k), and the fitting polynomial can be obtained:

[0105]

[0106] Where:

[0107]

[0108]

[0109]

[0110] P(i',j',k',) is the fitting value;

[0111] (x i' ,y j' ,z k' ) is the position of the element D(i',j',k') in the smoothing window;

[0112] a svt To fit the polynomial coefficients, write it in matrix form: P = AB, where P is the matrix composed of each element in the sliding window after fitting, A is the matrix composed of the least squares fitting polynomial parameters, that is, the weight coefficient matrix, and matrix B is the elements in the sliding window according to the coefficient a in the fitting polynomial. svt The corresponding x i' s y j' t z k' v A matrix of values;

[0113] l, m, and n represent the number of elements in the sliding window in the three directions of x, y, and z respectively; x i' s y j' t z k' v Represents the position of each point in the sliding window (x i' ,y j' ,z k' ) is the value calculated according to the polynomial fitting formula.

[0114] S32. The mean square error E=PD in the sliding window obtained by fitting can be determined, where the matrix E is the mean square error and the matrix D is a matrix composed of the initial values ​​of the elements in the sliding window;

[0115] S33. According to the least squares principle, the root mean square error is minimized, and the weight coefficient matrix A = B can be obtained by transforming the mean square error E. -1 D.

[0116] S4. The S determined in step S3 is X ×S Y ×S Z The weight matrix A is applied to S3 X ×S Y ×S Z All elements of the sliding window can obtain the best fitting value at the center point D(i,j,k) of the sliding window;

[0117] S5. Smoothing the entire initial three-dimensional velocity model along the X, Y, and Z directions of the initial three-dimensional velocity model respectively. After the smoothing process is completed, a smoothed velocity model is obtained; and the smoothed velocity model is output for migration processing.

[0118] The method of smoothing the entire three-dimensional velocity data volume along the three directions of X, Y and Z respectively in step S5 comprises the following steps:

[0119] Select the sliding window from the initial position and select the sliding window S in increments of 1 unit. X ×S Y ×S Z Time window, judge sliding S X ×S Y ×S Z Whether the end of the time window reaches the boundary of the initial three-dimensional velocity model in this direction, if not, repeat the operations of steps S3 and S4; in the X, Y, and Z directions of the initial three-dimensional velocity model, when the sliding S X ×S Y ×S Z When the end of the time window reaches the maximum element of the initial three-dimensional velocity model in this direction, the sliding is ended, and the least square convolution fitting weighted processing of the initial three-dimensional velocity model in this direction is stopped;

[0120] When processing in all three directions stops, the smoothing process of the entire initial three-dimensional velocity data volume is completed.

[0121] Example 1

[0122] The present invention discloses a velocity model smoothing processing method based on least square convolution fitting, which is used to process the overlying shale interbed target in T2 of a horizontal well area in Daqing Oilfield, and comprises the following steps:

[0123] S1. Combination select Select the migration method based on the geological requirements of the target block;

[0124] S2. Based on the selected migration method, the optimal sliding window parameter S relative to the target layer is calculated through the initial 3D velocity model. X , S Y , S Z ; The optimal sliding window parameter S is obtained X , S Y , S Z , determine S X ×S Y ×S Z Sliding window;

[0125] Where: S X is the smoothing window parameter of the main survey line direction of the 3D velocity model;

[0126] S Y is the smoothing window parameter of the contact line direction of the 3D velocity model;

[0127] S Z is the vertical smoothing window parameter of the 3D velocity model.

[0128] S21. This work area needs to solve the problem of optimizing the well structure error and clarify the underground depth of the geological structure layer, so the prestack depth migration PSDM processing method is selected. The prestack depth migration PSDM processing method is used to smooth the depth domain layer velocity model. First, the imaging target layer depth H is determined in meters, and the main line direction smoothing parameter S of the three-dimensional velocity model is calculated. X =4H / Δx+1, the smoothing window number S of the contact line direction of the three-dimensional velocity model Y =4H / Δy+1, vertical smoothing parameter S of the three-dimensional velocity model Z According to the demand for vertical smoothness of the model, the vertical smoothing parameter S is recommended. Z Take S X and S Y The larger value of S, where Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in units of m, S X , S Y , S Z Must be odd numbers, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness. Figure 1In the technical process, the processing object is the velocity model used for depth migration PSDM. The bin size in the observation system is 20m×10m, Δx=20, Δy=10, and the lower limit of the depth of the target layer is 2000m, then H=1000m. The main line direction smoothing parameter S of the 3D velocity model is X =4×1000 / 20+1=201, the smoothing window number S of the contact line direction of the three-dimensional velocity model Y =4×1000 / 10+1=401, S Z The possible value is 401.

[0129] S22. Combining geological requirements and migration method selection, the main problem in this work area is high-fidelity and high-precision imaging of the overlying mudstone interlayers of T2, so prestack time migration processing is selected.

[0130] The PSTM method of prestack time migration is used to smooth the RMS migration velocity model in the time domain. First, the imaging target layer time T (in seconds) and the average velocity V (in m / s) of the target layer are determined, and the main line direction smoothing parameter S of the three-dimensional velocity model is calculated. X =2VT / Δx+1, the smoothing window number S of the connection line direction of the three-dimensional velocity model Y =2VT / Δy+1, vertical smoothing parameter S of the three-dimensional velocity model Z According to the demand for vertical smoothness of the model, the vertical smoothing parameter S is recommended. Z Take S X and S Y The larger value of Δx and Δy is the size of the bin of the observation system of the processing object in the x and y directions. This parameter is included in the field shot point, detection point, and shot-check point relationship file formed during the seismic data acquisition process. Its meaning is the detection point spacing and shot point spacing. The unit is m, S X , S Y , S Z Must be odd numbers, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness. Figure 1 In the technical process, the processing object is the layer velocity model used for time migration PSTM. The bin size in the observation system is 20m×10m, Δx=20, Δy=10, the lower limit of the depth of the target layer is 1s, and the corresponding average formation velocity is 2000m / s. The main line direction smoothing parameter S of the 3D velocity model is X =2×1×2000 / 20+1=201, the smoothing window number S of the contact line direction of the three-dimensional velocity model Y=2×1×2000 / 10+1=401, S Z The value can be 401. S3. Substitute the S calculated in step S2 into X , S Y , S Z As the input sliding window parameter, the S X ×S Y ×S Z The speed values ​​of the sample points are fitted with a multivariate n-order polynomial, and the fitting error of the polynomial is minimized by the least squares principle to determine the S corresponding to the polynomial. X ×S Y ×S Z Weight coefficient matrix A, requiring polynomial degree p X , S Y , S Z , n is usually 3. The sliding window size calculated by S2 is S X ×S Y ×S Z (401x201x401), the number of elements in the x, y, and z directions are 401, 201, and 401 respectively, and the element at the center is D(201,101,201), and the fitting polynomial is:

[0131]

[0132] Where:

[0133] i'=i+x i' ;x i' =-200,-199,…,-1,0,1,…,199,200

[0134] j'=j+y j' ;y j' =-100,-99,…,-1,0,1,…,99,100

[0135] k'=k+z k' ;z k' =-200,-199,…,-1,0,1,…,199,200

[0136] P(i',j',k,) is the fitted value;

[0137] The fitting polynomial is written in matrix form, namely P=AB, and the weight coefficient matrix A is solved by calculating the minimum mean square error of the fitting value.

[0138] S4. The S determined in step S3 is X ×S Y ×S Z ​The weight coefficient matrix A is applied back to each element of the matrix composed of 401x201x401 elements of the initial velocity model, and the best fitting value P(201,101,201) at the center point D(201,101,201) of the sliding window can be obtained; for different sliding windows, there is a best fitting value P(i,j,k) corresponding to the center point value D(i,j,k) of its sliding window, and quality control can be performed by comparing the equally spaced sliding window fitting values ​​for the target layer with the actual values ​​of the corresponding points. The interval here is determined according to the velocity change trend of the target layer in the actual target block. For unconventional oil and gas reservoirs, 200 meters is selected. Through actual data testing in multiple unconventional oil and gas exploration areas in Daqing Oilfield, the sliding smoothing window parameters and the number of fitting polynomials selected by the method of the present invention are the optimal solutions.

[0139] S5. Slide along the X, Y, and Z directions in increments of 1 unit. X ×S Y ×S Z time window, repeat the operations of steps S3 and S4, in the three directions of the three-dimensional velocity model, when the sliding S X ×S Y ×S Z When the end of the time window reaches the maximum element of the data volume in that direction, the sliding stops, and the least square convolution fitting weighted processing of the velocity model in that direction stops. When the processing stops in all three directions, the smoothing processing of the entire three-dimensional velocity model ends. The smoothed velocity model is output for migration processing. Figure 2 It is the RMS velocity model profile after smoothing using the default parameters of the commercial seismic data processing software Omega. Figure 3 is the RMS velocity model profile after smoothing by the method of the present invention, Figure 2 Same location. Figure 4 is with Figure 2 , Figure 3 The actual seismic profile corresponding to the same position of the velocity profile. It can be seen that Figure 2 The smoothed velocity profile has multiple lateral velocity mutations and reversals at the 1500ms and 2000ms scale lines. Figure 2 The velocity profile after smoothing is Figure 4 There is no good correspondence between the actual underground structures, and the rapid lateral changes and obvious velocity boundaries will lead to false images and anomalies in the migration imaging. Figure 3 After smoothing by the method of the present invention, the velocity model profile as a whole does not show lateral velocity anomaly and reversal, and the velocity trend and Figure 4The actual underground structure is relatively consistent. A stable and correct velocity model is the basis of migration imaging. Figure 5(a) shows the calibration of the velocity model migration result section and the well synthetic record obtained by the conventional processing commercial software median filtering velocity model smoothing method. Figure 5(b) shows the calibration of the velocity model result section and the well synthetic record after the velocity model is smoothed by the method of the present invention. The correlation coefficient of the conventional velocity model smoothed migration section to the well synthetic record is about 0.75, and the wave resistance strength relationship is not good for the well. The correlation coefficient of the velocity model migration section obtained by the patent of the present invention to the well synthetic record is increased to 0.91, and the wave resistance strength relationship is significantly improved.

[0140] Seismic data processors and interpreters may make simple changes of one kind or another, or adjust and optimize or roughen the algorithm, or omit one or some steps in the process of the present invention according to the method disclosed in the present invention. Any practice that adopts the same core concept as the present invention shall be deemed as an infringement of this patent.

[0141] Note: Under the guidance of the specification published by the patent of this invention, researchers in various fields will take the core guiding ideas and research ideas published in this patent as guidance, refer to the adoption of all or part of the process, and make this or that change to obtain the same or similar results as all or part of the effects of the examples described in the specification, which should also be regarded as infringement.

Claims

1. A velocity model smoothing method based on least squares convolution fitting, characterized in that: The following steps are involved: S1. Select the target block and choose the migration method based on the geological requirements of the target block; S2. Based on the selected migration method, the optimal sliding window parameter S relative to the target layer is calculated through the initial 3D velocity model. X , S Y , S Z ; The optimal sliding window parameter S is obtained X , S Y , S Z , determine S X ×S Y ×S Z Sliding window; Where: S X is the smoothing window parameter of the main survey line direction of the 3D velocity model; S Y is the smoothing window parameter of the contact line direction of the 3D velocity model; S Z is the vertical smoothing window parameter of the 3D velocity model; S3. S calculated in step S2 X , S Y , S Z , for S in the sliding window X ×S Y ×S Z The velocity values ​​of the sample points are fitted with a multivariate n-order polynomial to determine the S corresponding to the polynomial. X ×S Y ×S Z Weight coefficient matrix A; S4. The S determined in step S3 is X ×S Y ×S Z The weight matrix A is applied to S3 X ×S Y ×S Z For all elements of the sliding window, obtain the best fitting value at the center point D (i, j, k) of the sliding window; S5. Smoothing the entire initial three-dimensional velocity model along the X, Y, and Z directions of the initial three-dimensional velocity model, after which a smoothed velocity model is obtained; outputting the smoothed velocity model for migration processing; The migration method includes: using a pre-stack depth migration processing method to smooth the depth domain layer velocity model and using a pre-stack time migration processing method to smooth the time domain migration velocity model; The method for smoothing the depth domain layer velocity model by using the prestack depth migration processing method comprises: First, determine the depth H of the imaging target layer; Calculate the smoothing window parameter S of the main survey line direction of the 3D velocity model X ; And the smoothing window parameter S of the contact line direction of the 3D velocity model Y ; and the vertical smoothing window parameter S of the 3D velocity model Z ; Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness; and / or, S X =4H / Δx+1; S Y =4H / Δy+1; S X , S Y , S Z All are odd numbers; Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in meters; H is the depth of the imaging target layer, in meters; The vertical smoothing window parameter S of the three-dimensional velocity model Z , the vertical smoothing window parameter S is determined according to the requirements of the vertical smoothness of the model. Z Take S X and S Y The larger value among The method of smoothing the time domain migration velocity model by using the prestack time migration processing method comprises: First, determine the imaging target layer time T and the average velocity V of the target layer; Calculate the smoothing window parameter S of the main survey line direction of the 3D velocity model X ; And the smoothing window parameter S of the contact line direction of the 3D velocity model Y ; and the vertical smoothing window parameter S of the 3D velocity model Z ; Relative to the elements in the sliding window of the input, S X , S Y , S Z The larger the value, the higher the smoothness. X , S Y , S Z The smaller the value, the lower the smoothness; and / or, S X =2VT / Δx+1; S Y =2VT / Δy+1; S X , S Y , S Z All are odd numbers; Where: Δx and Δy are the sizes of the observation system surface element of the processing object in the x and y directions, in units of m; T is the imaging target layer time, in seconds; V is the average velocity of the target layer, in m / s; The S Z The vertical smoothing window parameter S is determined according to the requirement for vertical smoothness of the model. Z Take S X and S Y The larger value.

2. The velocity model smoothing processing method based on least squares convolution fitting according to claim 1 is characterized in that: In step S3, X ×S Y ×S Z The method for obtaining the weight coefficient matrix A is: S71. Determine the least squares fitting polynomial; S72. Determine the mean square error within the sliding window obtained by fitting; S73. According to the least squares principle, the root mean square error is minimized, and the weight coefficient matrix A can be obtained by transforming the mean square error E.

3. The velocity model smoothing processing method based on least squares convolution fitting according to claim 2 is characterized in that: In step S71, the method for determining the least square fitting polynomial is: S X , S Y , S Z The number of elements in the x, y, and z directions in the sliding window of size is l, m, and n respectively, and the element at the center is D (i, j, k). The fitting polynomial is: ； Where: ； is the fitted value; To smooth the elements in the window Position within the smoothing window; a svt To fit the polynomial coefficients, write it in matrix form: P = AB, where P is the matrix composed of each element in the sliding window after fitting, A is the matrix composed of the least squares fitting polynomial parameters, that is, the weight coefficient matrix, and matrix B is the elements in the sliding window according to the coefficient a in the fitting polynomial. svt Corresponding A matrix of values; l, m, and n represent the number of elements in the three directions of xyz contained in the sliding window; Represents the position of each point in the sliding window The value calculated according to the polynomial fitting formula; and / or, The step S72 determines the mean square error E=PD in the sliding window obtained by fitting; Among them: Matrix E is the mean square error; Matrix D is the matrix composed of the initial values ​​of each element in the sliding window; and / or, S73 obtains the weight coefficient matrix A=B by transforming the mean square error E. -1 D.

4. The velocity model smoothing processing method based on least squares convolution fitting according to claim 1 is characterized in that: In step S3, the polynomial degree n X 、n Y 、n Z , n is 3.​​​ 5. The velocity model smoothing processing method based on least squares convolution fitting according to claim 1 is characterized in that: The method of step S5 for smoothing the entire initial three-dimensional velocity model along the three directions of X, Y, and Z of the initial three-dimensional velocity model respectively comprises: Select the sliding window from the initial position and select the sliding window S in increments of 1 unit. X ×S Y ×S Z Time window, judge sliding S X ×S Y ×S Z Whether the end of the time window reaches the boundary of the initial three-dimensional velocity model in this direction, if not, repeat the operations of steps S3 and S4; in the X, Y, and Z directions of the initial three-dimensional velocity model, when the sliding S X ×S Y ×S Z When the end of the time window reaches the boundary of the initial three-dimensional velocity model in that direction, the sliding is terminated, and the least square convolution fitting weighted processing of the initial three-dimensional velocity model in that direction is stopped; When processing in the X, Y, and Z directions stops, the smoothing process of the entire initial three-dimensional velocity model is completed.

Images