Near-surface longitudinal wave velocity conversion method based on static sounding data

Abstract

The invention discloses a near-surface longitudinal wave velocity conversion method based on static sounding data, and the method comprises the steps: building a nonlinear mapping relation among conical tip resistance, side friction resistance and longitudinal wave velocity through an artificial neural network; and the advantages of convenient layout and low cost of the static sounding method are further combined, so that the data volume and the acquisition efficiency of site longitudinal wave velocity survey are considered, and the construction cost is reduced. The application steps of the method comprise four stages of work area training set creation and preprocessing, initial artificial neural network establishment, artificial neural network training and artificial neural network generalization. A site blind well test experiment proves that the shallow rock-soil body longitudinal wave velocity conversion precision of the static sounding-longitudinal wave velocity conversion method based on the artificial neural network reaches 90% or above, and compared with a traditional empirical formula method, the static sounding-longitudinal wave velocity conversion method based on the artificial neural network has higher accuracy, longitudinal resolution and robustness, and is beneficial to discrimination of a low speed reduction zone and a ghost reflection interface.

Description

Technical Field

[0001] The present invention belongs to the field of rock physics research, and involves artificial neural network technology, data preprocessing technology and global optimization search algorithm. The invention results serve to provide an important basis for determining near-surface lithology, dividing the surface low-velocity reduction zone, and determining the virtual reflection interface, as well as providing theoretical guidance for the optimization of stimulating lithology and well depth, and providing basic parameters and geological understanding for shallow correction in data processing, so as to improve the quality of seismic exploration data. Background Art

[0002] The spatial longitudinal and lateral variation of near-surface P-wave velocity is the key target of near-surface structure investigation. It is an important basis for identifying near-surface lithology, dividing the surface low-velocity reduction zone, and determining the virtual reflection interface. It can provide theoretical guidance for the optimization of stimulating lithology and well depth, and provide basic parameters and geological understanding for shallow correction in data processing. Therefore, the development of a fast and efficient P-wave velocity structure investigation method is the key foundation for improving the quality of exploration data and reducing multi-solution.

[0003] However, conventional direct velocity testing methods (such as drilling coring, micro-well logging, small refraction, and cannon first-arrival tomography inversion) are often limited by construction costs and construction cycles, and it is difficult to achieve large-scale coverage in practice. Static penetration is an in-situ testing technology that uses static force to press the probe into the soil layer to obtain the physical and mechanical properties of the soil. The advantages of this method are: it is not limited by the wavelength of elastic waves and the effective frequency band of seismic detectors, and has the advantages of high resolution, simple construction, strong immediacy, and low cost. It is widely used in engineering for rock and soil stratification, judging sand liquefaction, bedrock weathering, and other purposes. The limitation of this method is that the readings are recorded as cone tip resistance and side friction resistance information, which reflects the distribution of compression and shear modulus of the rock and soil body, but cannot directly measure velocity information.

[0004] Therefore, if the static penetration test and velocity response can be integrated to serve the site wave velocity structure analysis together, it is possible to take into account the test efficiency and reduce the construction cost on the basis of the accuracy of conventional survey data. To this end, some scholars have tried to establish a linear formula between the cone tip resistance and the compression modulus through statistical laws to associate the static penetration test response with the soil wave velocity, and have achieved certain application effects in practice. However, the limitations of this method are also significant. First, the existing empirical formulas are mostly based on linear or exponential relationship assumptions, so the structure and physical state of the soil have to be simplified, so they cannot reflect the complex rock physics relationship and the scope of application is limited; second, the fitting coefficient is highly regional and cannot be adjusted with the change of the measurement point location and depth; third, some models based on wave theory involve too many parameters to be obtained, and a large number of indoor experiments are required, which cannot achieve the purpose of cost saving. Summary of the invention

[0005] The present invention aims to address the problems of high cost and long cycle of traditional direct measurement methods, as well as the applicability difficulties faced by the existing empirical formula method for wave velocity conversion. A multi-source data fusion algorithm based on artificial neural networks is developed. Through a large number of sample supervised learning, the nonlinear mapping relationship between cone tip resistance, side friction resistance, depth and wave velocity is learned to take into account the accuracy of micro-logging data and the coverage area of ​​static penetration data. A corresponding implementation process is proposed for investigation practice, in order to form a high-efficiency, low-cost, and high-precision near-surface fine velocity structure investigation model.

[0006] Technical solution:

[0007] A near-surface P-wave velocity conversion method based on static penetration data comprises the following steps:

[0008] S1. Creation and preprocessing of the training set of the work area, where the data set is composed of micro-logging interpretation speed v p The static penetration data are constructed in pairs according to the depth of the measuring points;

[0009] S2, initial artificial neural network establishment;

[0010] S3, artificial neural network training;

[0011] S4. Obtain the static penetration point to be converted, and input the preprocessed point into the trained artificial neural network model to obtain the velocity variation curve of the point with depth.

[0012] Preferably, the static penetration test data includes the cone tip resistance q c and lateral friction f s .

[0013] Preferably, the architecture of the initial artificial neural network in S2 is: combined cone tip resistance q c , lateral friction resistance f s With depth d, construct three-dimensional feature space, longitudinal wave velocity v p as target output.

[0014] Preferably, the initial artificial neural network in S2 adopts a multi-layer feedforward learning mechanism, and its matrix expression is as follows:

[0015]

[0016] Where h(x)=1 / (1+e -x ) indicates that the result of weighted summation of neurons in each layer is passed to the next layer through the sigmoid activation function; l indicates the number of hidden layers; W (k) It represents the matrix form of the weights connecting the kth layer in the hidden layer with all neurons in the previous layer. The connection weight matrix between the output layer and the last layer of the hidden layer is recorded as W(o) ; Basic threshold B (k) is a row vector whose number of elements is equal to the number of neurons in this layer, where the element is denoted by b i (k) ; A (k) represents the output of the k-th layer of neurons; σ is the activation function of the output layer; is the output of the artificial neural network.

[0017] Preferably, the purpose of training in S3 is to obtain a set of network connection weights W and basic thresholds B so that the neural network output Compared with the actual observation result Y n The mean square error is the smallest.

[0018] Preferably, the specific steps of training in S3 are:

[0019] ① Given the initial weight matrix and the initial basic threshold B, the initial model x is generated k =x 0 , set the target error ε>0, set k=0, give the step size enlargement factor α>1 and the step size reduction factor 0<β<1;

[0020] ② It is stipulated that when a set of data is input, the corresponding neural network outputs once Compared with the actual observation result Y n The difference satisfies the functional relationship f(x), that is: Then we can calculate the Jacobian matrix of f(x), that is, Guide the update of model parameters;

[0021] ③ Obtain the model update amount Δx k+1 ;

[0022] ④ Update x according to the following formula k+1 , x k+1 =x k -Δx k+1 ;

[0023] ⑤If S(x k+1 )≥S(x k ), then let λ = αλ, and return to step ② to loop, where S is the objective function in the sense of least squares, that is, the overall performance of the model in all training:

[0024] ⑥ Update x k , let λ=βλ; if ||A(x k ) T f(x k )|| 2 ≤ε, the iteration terminates, where ε is a small positive number; otherwise, it returns to step ② and continues the loop.

[0025] Preferably, the model update amount in step ③ is obtained by the following formula:

[0026] Δx k+1 =(A(x k ) T A(x k )+λ k I) -1 A(x k ) T f(x k )

[0027] Where I represents the identity matrix and the coefficient λ acts as follows: k ) T A(x k )) -1 A(x k ) T f(x k )|| 2 ≤ε, with λ k =0; otherwise λ k >0, adjust λ k So that A(x k ) T A(x k )+λ k I is positive.

[0028] Beneficial effects of the present invention

[0029] The velocity conversion method proposed in the present invention can obtain effective shallow longitudinal wave velocity structure from static penetration data when micro-well logging is difficult to implement, the deployment density is low, or the construction period is limited, which is conducive to reducing the cost of field investigation. In principle, this technology does not rely on theoretical assumptions and indoor tests, so it has good applicability to sites with different geotechnical properties. Comparative tests show that the static penetration data conversion method based on artificial neural network has a conversion accuracy of more than 90% for shallow velocity, which has higher accuracy, longitudinal resolution and robustness than the traditional empirical formula conversion method. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 This is a network diagram for longitudinal wave velocity conversion using a multi-layer feedforward artificial neural network model.

[0031] Figure 2 This is a flow chart of the near-surface P-wave velocity conversion algorithm based on data fusion technology.

[0032] Figure 3 This is the distribution map of the near-surface velocity survey test points.

[0033] Figure 4This is the longitudinal wave velocity map obtained based on the artificial neural network model and empirical formula method. DETAILED DESCRIPTION

[0034] The present invention will be further described below in conjunction with embodiments, but the protection scope of the present invention is not limited thereto:

[0035] Conventional direct velocity measurement methods are often limited by site conditions and construction costs, and it is difficult to ensure both data volume and data accuracy. The present invention takes advantage of the convenience and low cost of static penetration testing and proposes a method of fusing formation resistance and velocity information through an artificial neural network to achieve reliable shallow longitudinal wave velocity structure acquisition. The actual process of this technology includes three steps: data preparation, model training, and model generalization: based on the cone tip resistance, side friction resistance, and deep structural feature space, a nonlinear mapping relationship between data is established through a multi-layer feedforward training mechanism to obtain the optimal model, and finally the conversion velocity curve is output.

[0036] Combination Figure 2 , which specifically includes the following implementation steps:

[0037] (1) Creation and preprocessing of work area training set:

[0038] Traditional micro-logging interpretation often only focuses on layer velocity, while ignoring the velocity changes between layers. Therefore, there is a mismatch between the vertical resolution of the micro-logging and the two cannot be directly paired for training. Therefore, the present invention intends to extract the velocity information of the formation based on the first arrival time of the micro-logging, that is, to arrange the arrival time of each shot from top to bottom according to the depth, obtain the arrival time record of the point and calculate the slope as the velocity value of the point, and then linearly smooth the 5 consecutive sampling points in the depth domain to reduce the interference of local data singular points on the overall training effect. Subsequently, the nearest micro-logging interpretation velocity (v p ) and static penetration data (cone tip resistance q c and lateral friction f s ) Pair the two points according to the depth of the measuring points. During this period, the apparent depth must be projected and corrected according to the starting elevation and field information such as well inclination to ensure that the paired data are from the same horizontal plane. That is, the interpretation, calibration and pairing of micro-logging and static penetration data are completed.

[0039] (2) Initial artificial neural network establishment:

[0040] Assuming the number of paired data is n, the combined cone tip resistance q c , lateral friction resistance f s With depth d, construct three-dimensional feature space, longitudinal wave velocity v p As the target output, it is recorded as:

[0041]

[0042] Combination Figure 1 , using a multi-layer feedforward learning mechanism as the training logic, the parameters are set as follows: the number of input parameters is 3, the number of output layers is 1; the number of hidden layers is l, and the number of neurons in each layer is m. The matrix form of the weights of the kth layer in the hidden layer connected to all neurons in the previous layer (the first layer of the hidden layer is connected to the input layer) is W (k) The number of rows in the matrix is ​​the same as the number of neurons in this layer, the number of columns is the same as the number of neurons in the previous layer, and the element is the weighted value w connecting the i-th neuron with the j-th neuron in the previous layer. ij (k) ; Basic threshold B (k) is a row vector whose number of elements is equal to the number of neurons in this layer, where the element is denoted by b i (k) ; The result of weighted summation of neurons in each layer is passed to the next layer through the sigmoid activation function, and its expression is:

[0043] h (x)＝ 1 / (1 +e -x ) (2)

[0044] At the same time, the connection weight matrix between the output layer and the last layer of the hidden layer is recorded as W (o) , where the weighted value of the connection between the element i-th neuron and the j-th neuron in the last layer of the hidden layer is w ij (o) The output required by the present invention has only speed attribute, namely W (o) Degenerates into a column vector, the number of elements is the number of neurons in the last hidden layer m. Basic threshold B (o) The element is denoted by b i (o) Based on this, the multi-layer feedforward artificial neural network constructed by the present invention can be expressed in the following matrix form:

[0045]

[0046] Wherein σ is the activation function of the output layer, and the present invention adopts the identity function, that is, σ(x)=0.

[0047] (3) Artificial neural network training:

[0048] The purpose of training is to find a set of network connection weights W and basic thresholds B so that the neural network output Compared with the actual observation result Y n The mean square error is the smallest. The mean square error function E is defined as:

[0049]

[0050] The feature space is randomly divided into 10 equal parts and used as validation sets in turn. The target optimization function is to minimize the average prediction error of 10 validations. The connection weight matrix and the basic threshold matrix are obtained by searching through the back propagation algorithm. The training process is: the number of weights from the input layer to the hidden layer is a×m1, and the number of weights between hidden layers is The number of weights from the hidden layer to the output layer is b×m l , the total number of thresholds is The sum of the number of hidden layer neurons and the number of output layer neurons b. Then, if the unknown number is represented by the vector x = [x1 x2 … x M ] indicates that:

[0051]

[0052] Rewrite the objective function as:

[0053] minS(x)=f T (x)f(x)=||f(x)|| 2 , (5)

[0054] Assume that there are N sets of training data, where f(x) = [f1(x) f2(x) … f N (x)], the inverse problem can be transformed into a nonlinear least squares problem and solved by continuous iteration. Taking the Levenberg-Marquardt algorithm as an example, f(x) is replaced by x at the kth iteration point. k Perform Taylor first-order expansion and specify x k The trust region of is ε, that is:

[0055]

[0056] Then the solution to this equation satisfies:

[0057] (A(x k ) T A(x k )+λ k I)z=-A(x k ) T f(x k ) (7)

[0058] So as to satisfy the recursive relationship,

[0059] x k+1 =x k -(A(x k ) T A(x k )+λ k I) -1 A(xk ) T f(x k ) (8)

[0060] When ||(A(x k ) T A(x k )+λ k I) -1 A(x k ) T f(x k )|| 2 ≤ε, with λ k =0; otherwise λ k >0, adjust λ k Make

[0061] A(x k ) T A(x k )+λ k I is positive definite, so that the z generated by equation (8) is in the descending direction.

[0062] In summary, the model parameter iteration process is as follows:

[0063] ① Given the initial weight matrix and threshold vector B, the initial model x is generated k =x 0 , set the target error ε>0, set k=0, give the step size magnification factor α>1 and the step size reduction factor 0<β<1; ② calculate the Jacobian matrix of f(x) ② Calculate the Jacobian matrix of f(x) ③ Solve equation (8) to get the model update; ④ Update x according to the following formula k+1 , x k+1 =x k -Δx k+1 ⑤If S(x k+1 )≥S(x k ), then set λ = αλ and return to step ② to loop; ⑥ Update x k , let λ=βλ; if ||A(x k ) T f(x k )|| 2 ≤ε, the iteration terminates, otherwise returns to step ② to continue the loop.

[0064] (4) Generalization of artificial neural networks:

[0065] The trained artificial neural network model is applied to unpaired static penetration points to predict the velocity change curve of a single well and the velocity structure of the work area. First, static penetration points that have not participated in the pairing training are selected in the work area. After preprocessing, they are input into the trained artificial neural network model to obtain the velocity change curve of the point with depth. The conversion effect is quantitatively evaluated by the correlation coefficient R:

[0066]

[0067] The superscript * indicates the micro-logging extraction speed, and the closer R is to 1, the better the conversion effect.

[0068] like Figure 3 , taking a certain work area as an example, the study was carried out, and the speed conversion test of the conversion effect was carried out through paired static penetration-micro logging points, and the effectiveness of the algorithm of the present invention was verified by comparing it with the traditional conversion method. A total of 4 test points were set up in the experiment. Each point was simultaneously subjected to a double-well micro logging survey and a double-bridge static penetration test. The micro logging depth was 30m, and the static penetration depths from S1 to S4 were 30.3m, 27.2m, 24.3m and 30.3m respectively. The near-surface lithology was mainly clay, silty clay and fine sand. The output results before and after dimensionality reduction were compared by blind well verification, that is, only three point data (S1, S2, S4) were used as training samples, and the remaining point (S3) was used as a verification point and did not participate in the training process. The initial neural network model used a structure with 10 hidden layers and 10 neurons in each layer.

[0069] The method of the present invention is compared with the traditional static penetration test-longitudinal wave velocity empirical formula method. Taking the common exponential conversion relationship as an example, this method believes that the longitudinal wave velocity satisfies:

[0070]

[0071] The same blind well test method is used. Substituting the above formula into S1, S2, S4, the result is as follows Figure 4As shown by the gray curve in . Obviously, this conversion formula can only describe the increasing trend of velocity with depth, and is not as good as the conversion accuracy of the method of the present invention, which is specifically manifested as follows: ① After calculation, the regression coefficient R between the converted P-wave velocity (gray curve in the figure) and the micro-logging interpretation velocity (red curve in the figure) obtained by the traditional empirical formula method is only 0.78, and the longitudinal resolution is poor, which can only reflect the general growth trend, but the ability to describe the local detail changes is insufficient; ② The robustness of the traditional empirical formula method is poor, that is, it is too sensitive to local outliers in static penetration, resulting in a large-scale data drift problem; ③ After calculation, the regression coefficient R between the result (blue curve) obtained by the velocity conversion method proposed by the present invention and the micro-logging interpretation velocity (red curve in the figure) reaches 0.91, which is significantly better than the empirical formula method, and has higher resolution in the depth domain, stronger robustness to outliers, and no data drift phenomenon occurs.

[0072] The specific embodiments described herein are merely examples of the spirit of the present invention. Those skilled in the art may make various modifications or additions to the specific embodiments described or replace them in similar ways, but they will not deviate from the spirit of the present invention or exceed the scope defined by the appended claims.

Claims

1. A near-surface longitudinal wave velocity conversion method based on static penetration data, characterized in that It includes the following steps: S1. Creation and preprocessing of the training set of the work area, where the data set is composed of micro-logging interpretation speed v p The static penetration data are constructed in pairs according to the depth of the measuring points; S2, initial artificial neural network establishment; S3, artificial neural network training; S4. Obtain the static penetration point to be converted, and input the preprocessed point into the trained artificial neural network model to obtain the velocity variation curve of the point with depth.

2. The method according to claim 1, characterized in that The static penetration test data includes the cone tip resistance q c and lateral friction f s .

3. The method according to claim 2, characterized in that The architecture of the initial artificial neural network in S2 is: combined cone tip resistance q c , lateral friction resistance f s With depth d, construct three-dimensional feature space, longitudinal wave velocity v p as target output.

4. The method according to claim 1, characterized in that: The initial artificial neural network in S2 adopts a multi-layer feedforward learning mechanism, and its matrix is ​​expressed as follows: Where h(x)=1 / (1+e -x ) indicates that the result of weighted summation of neurons in each layer is passed to the next layer through the sigmoid activation function; l indicates the number of hidden layers; W (k) It represents the matrix form of the weights connecting the kth layer in the hidden layer with all neurons in the previous layer. The connection weight matrix between the output layer and the last layer of the hidden layer is recorded as W (o) ; Basic threshold B (k) is a row vector whose number of elements is equal to the number of neurons in this layer, where the element is denoted by b i (k) ; A (k) represents the output of the k-th layer of neurons; σ is the activation function of the output layer; is the output of the artificial neural network.

5. The method according to claim 4, characterized in that The purpose of training in S3 is to find a set of network connection weights W and basic thresholds B so that the neural network output Compared with the actual observation result Y n The mean square error is the smallest.

6. The method according to claim 4, characterized in that The specific steps of training in S3 are: ① Given the initial weight matrix and the initial basic threshold B, the initial model x is generated k =x 0 , set the target error ε>0, set k=0, give the step size enlargement factor α>1 and the step size reduction factor 0<β<1; ② It is stipulated that when a set of data is input, the corresponding neural network outputs once Compared with the actual observation result Y n The difference satisfies the functional relationship f(x), that is: Then we can calculate the Jacobian matrix of f(x), that is, Guide the update of model parameters; ③ Obtain the model update amount Δx k+1 ; ④ Update x according to the following formula k+1 , x k+1 =x k -Δx k+1 ; ⑤If S(x k+1 )≥S(x k ), then let λ=αλ, and return to step ② to loop, where S is the objective function in the sense of least squares; ⑥ Update x k , let λ=βλ; if ||A(x k ) T f(x k )|| 2 ≤ε, the iteration terminates, where ε is a small positive number; otherwise, it returns to step ② and continues the loop.

7. The method according to claim 6, characterized in that The model update amount in step ③ is obtained by the following formula: Δx k+1 =(A(x k ) T A(x k )+λ k I) -1 A(x k ) T f(x k ) Where I represents the identity matrix and the coefficient λ acts as follows: k ) T A(x k )) -1 A(x k ) T f(x k )|| 2 ≤ε, with λ k =0; otherwise λ k >0, adjust λ k So that A(x k ) T A(x k )+λ k I is positive.

Images