The complex bedding structure of shale may make its ultrasonic characteristics complicated, and make it difficult to calculate the mechanical parameters of the shale reservoir accurately. Analyzing the data of shale gas blocks, Josh et al.^[1] found that the complex bedding structure of the shale led to anisotropy of the acoustic logging data. Vernik et al.^[2] and Deng et al.^[3] reached similar conclusion through the anisotropy of shale velocity. Horne et al.^[4] used dipole acoustic logging data to judge the elastic anisotropy of layered shale, and the estimated anisotropy was consistent with the predicted value of adjacent wells. Zhang et al.^[5] measured the vertical and horizontal ultrasonic velocities of two rock masses, and found obvious velocity anisotropy parallel with and perpendicular to the bedding. Yan et al.^[6], Cheng et al.^[7,8] studied the wave velocity changes of shale under uniaxial and triaxial compression conditions respectively, and found that the P and S wave velocities increased first and then decreased with the increase of axial load. Chen et al.^[9], Xiong et al.^[10] and Wu et al.^[11] found through ultrasonic transmission experiments that the wave velocity decreased with the increase of the bedding angle or the decrease of the frequency, the attenuation coefficient increased with the increase of the angle or frequency, the velocity parallel with the bedding direction was positively correlated with the bedding density, while the velocity perpendicular to the bedding direction was negatively correlated with the bedding density. Schn et al.^[12] established a function of shale elastic properties based on ultrasonic experiments and calculated velocity and anisotropy by this function.

In summary, domestic and foreign researchers mainly used outcrops or downhole cores to make physical models with different bedding structures, and summarized the qualitative understanding of shale acoustic responses by ultrasonic transmission experiments^[13,14]. Since shale samples are difficult to prepare and a large number of shale samples are needed for and damaged in experimental study, it is impossible to carry out large-scale physical model experiments, so the research results have some limitations. In view of the above problems, taking the shale bedding structure as the research object, based on the wave theory, the ultrasonic transmission experiment process was numerically simulated to find out the ultrasonic response law of different layered models, and the gray system theory was adopted to screen the acoustic parameters sensitive to bedding structures, to establish a dynamic and mechanical shale model. The research results are of great significance for accurately predicting mechanical parameters with acoustic data.

The elasticity theory mainly studies the relationship between the applied force and deformation of an object. The wave equation is jointly established by the constitutive equation, motion equilibrium differential equation and geometric equation. For an elastic isotropic medium, the two-dimensional first-order velocity, stress elastic wave equation is:

In time, the approximate solution of the time difference equation with second-order precision is used, and in space, the approximate solution of the fourth-order difference equation is used. Using the staggered grid format, it is:

The difference format of the wave equation can be derived:

1.2.1. Boundary conditions

In the calculation area, the left and right boundaries are basically consistent with actual conditions, and are treated as conventional reflection boundaries, while absorption conditions are set at the upper and lower boundaries and the four corner points.

1.2.2. Initial conditions

When t≤0, the velocity and displacement are zero, i.e.

1.2.3. Source conditions

1.2.4. Convergence conditions

Analysis of error propagation and accumulation is the basis for ensuring the convergence of the difference format. Dong et al.^[15] deduced in detail the stability conditions of high-order difference equations for staggered grids. When the medium is uniform and isotropic, the time precision is second order, and the space precision is fourth order, the stability condition is:

Laboratory ultrasonic testing is measured using the transmission method. However, due to the heterogeneity of shale and the uncertainty of bedding distribution, it is difficult to study the influence of the bedding structure on the ultrasonic propagation characteristics. Moreover, unsmooth contact of the probe with the core and inaccurate initial time reading manually in the transmission experiment may cause errors too. Taking the vertical section of the measured core as the model, and the shot signal of the ultrasonic probe as the source, this numerical simulation derived the first-order displacement- stress elastic wave equations based on the wave theory, and calculated the half-length staggered grid and simulated the ultrasonic transmission experiments by the Matlab program.

1.3.1. Waveform comparison

In order to make the model consistent with the ultrasonic transmission experiment, the probe frequency was set at 250 kHz. The incident wave signal was displayed as an initial waveform. When the incident signal passed through the core and reached the receiving end, the laboratory obtained the transmitted signal, which was displayed as an experimental waveform, and the result of the numerical calculation was displayed as a simulated waveform (Fig. 1). The waveform comparison show that both the simulated waveform and the experimental waveform reached the receiving end at around 0.9×10^-5 s, the voltage for them were between -5 and 5 V in the similar range (Fig. 1).

Setting the simulated waveform signal as X, and the experimental waveform signal as Y, and the correlation of the two waveforms is calculated as:

The standard deviations of X and Y are:

The covariance is:

The correlation coefficient between the simulated waveform and the experimental waveform calculated is 0.85, indicating high coincidence between them. That means that the numerical simulation method can effectively simulate the ultrasonic transmission experiment.

1.3.2. Numerical calculation

Since the bedding angle is the most significant structural feature of shale, ideal core models with different bedding angles were established based on actual cores, and numerical simulation calculations were carried out on the models. Fig. 2 shows outcrop cores from the Lower Silurian Longmaxi Formation in southeastern Chongqing, coded A₁, A₂, A₃, and A₄. The bedding angle is the angle between the bedding direction and the propagating direction of ultrasonic wave. The bedding angle of No. A₁ A_2, A_3, and A₄ are 0°, 30°, 60°, and 90° respectively (Fig. 2). In order to further verify the validity of the numerical simulation method, numerical models B₁, B₂, B₃, B₄ were established by using the cores A₁, A₂, A₃ and A₄ as templates (Fig. 3). In the models, the velocity of the ultrasonic wave in the bedding was set at 1 900 m/s, and in the skeleton at 4 200 m/s. The models have ideal bedding angles, and correspond one-by- one to the cores A₁, A₂, A₃, and A₄. On that basis, the ultrasonic transmission experiment and numerical simulation results were compared to verify the numerical simulation method.

The attenuation coefficients of cores A₁, A₂, A₃ and A₄ at different probe frequencies were calculated by ultrasonic transmission experiments. The results show that at the same bedding angle, the attenuation coefficient increases as the probe frequency increases (Fig. 4a). At the same time, the attenuation coefficients of the models B₁, B₂, B₃ and B₄ at the same conditions were calculated by numerical simulation (Fig. 4b). The law obtained from the numerical simulation is basically consistent with the law from the ultrasonic transmission experiments^[16], that is, at the same bedding angle, the attenuation coefficient is positively correlated with the probe frequency.

The above numerical calculations are based on ideal bedding angle models. To further verify the reliability of the above numerical simulation results, numerical image processing method was used to extract the bedding angle of real cores, and then the attenuation coefficient at different bedding angles were calculated. Fig. 5 shows the photograph of outcrop cores from the Lower Silurian Longmaxi Formation in southeastern Chongqing, coded C₁, C₂, C₃, C₄. The bedding angle of C_1, C_2, C₃ and C₄ are 0°, 30°, 45°, and 90° respectively (Fig. 5).

Taking cores C₁, C₂, C₃ and C₄ as research objects, edge extraction and Hough line detection were carried out on the bedding and matrix skeleton in the photographs^[17], to establish four numerical simulation models, coded D₁, D₂, D₃, D₄ (Fig. 6) corresponding to C₁, C₂, C₃, C₄ in bedding angle, respectively. Since the digital image is composed of rectangular pixels, and each pixel is a small square, these small squares can be used as quadrilateral elements in finite difference calculation, so the digital image extracted for the bedding structure can be easily converted into finite difference grids by mapping the square elements. First, the coordinates of the four nodes of each unit were converted into the physical coordinates of the corresponding vector space according to the proportional relationship between the actual and pixel sizes of the image, and then the entire image could be converted into square finite difference grids. This method is called square unit mapping method. Then, according to the gray value of each pixel point, each unit was given a corresponding parameter, thereby the bedding structure characterized by the digital image was introduced to the numerical model. Finally, the Matlab algorithm was used to carry out numerical simulation calculation, and the calculation results were compared with the results of ultrasonic transmission experiments. The results show that the attenuation coefficient from the numerical simulation is positively correlated with the bedding angle, which is basically consistent with the law of the attenuation coefficient from the ultrasonic transmission experiment on cores (Fig. 7).

In summary, to cope with the high cost of ultrasonic transmission experiments, and errors caused by unsmooth contact between the core and the probe and human^[18], the results from numerical simulation method for shale ultrasonic transmission experiments proposed in this study are in good agreement with actual ultrasonic transmission experiment results. Therefore, the numerical simulation method can be used to study the acoustic characteristics of shale bedding structure.

A prominent feature of shale is developed bedding structure. The parameters characterizing bedding structure mainly include bedding angle, bedding thickness and bedding density. In order to carry out numerical simulation consistent with the conditions of ultrasonic transmission experiments, the core plugs in the models were set at 50 mm high and 25 mm in diameter, the P wave velocity was 4 200 m/s and the S wave velocity 2 500 m/s in the shale skeleton, the P wave velocity was 1 900 m/s and the S wave velocity 1 100 m/s in the shale bedding, the skeleton density was 2.6 g/cm³, and the bedding density was 2.3 beddings/cm². The incident wave from the laboratory probe was taken as the initial wave, and the direction of the ultrasonic wave kept unchanged. The cross section of the core was divided into grids (250 mm long and 125 mm wide) with the grid spacing of 0.2 mm and the time step of 20 ns. By changing the bedding model parameters, different bedding structures of shale were simulated to find out the influence of the bedding structure on ultrasonic characteristics.

2.1.1. Influence of bedding thickness on ultrasonic characteristics

When setting the model parameter, the bedding that is straight and parallel to the layer was defined as a horizontal bedding, and the bedding that is straight and perpendicular to the layer was defined as vertical bedding. At the same the bedding angle and density, the vertical bedding thickness was set at 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, and 2.0 mm respectively. The models of this group were coded E₁ to E₁₁. In another group of models, the horizontal bedding thickness was set at these thicknesses and coded F₁ to F₁₁ respectively. In the two groups of models, the ultrasonic transmitted from the bottom to the top of the core. In this study, the E₂, E₁₀ and F₂ and F₁₀ models with obvious bedding thickness and structure characteristics were selected to analyze the influence of bedding thickness on the ultrasonic characteristics, and analyze the effect of bedding thickness on ultrasonic characteristics microscopically further by wave field snapshot. The bedding thickness of E₂ and E₁₀ are 0.2 mm and 1.8 mm respectively, that of F2 and F₁₀ are 0.2 mm and 1.8 mm respectively too. From the numerical calculations of wave velocity, with the increase of bedding thickness, the wave velocity decreases linearly both in the horizontal bedding and the vertical bedding, but that in the horizontal bedding decreases more rapidly (Fig. 8a). The calculated results of attenuation coefficient show that with the increase of bedding thickness, the attenuation coefficient increases in both horizontal bedding and vertical bedding models, the attenuation coefficient in horizontal bedding increases in a power function, while that in the vertical bedding increases in an exponential function (Fig. 8b). The calculated results of dominant ultrasonic frequency and dominant amplitude show that, with the change of bedding thickness, the dominant ultrasonic frequency in the vertical bedding and horizontal bedding, and the dominant ultrasonic amplitude in the vertical bedding change irregularly, while the dominant ultrasonic amplitude in the horizontal bedding tends to decrease on the whole as the bedding thickness increases (Fig. 8c).

In the previous preset models, the vertical bedding thickness of E₁₀ was larger than the vertical bedding thickness of E₂, and the ultrasonic wave propagated from the bottom to the top of the core in both. The wave field snapshots show the distance that the ultrasonic signal can reach in E₁₀ is shorter than that in E₂ in the same time interval (Fig. 9a, 9b), so the ultrasonic wave velocity is lower in E₁₀ (Fig. 8a). The horizontal bedding propagation of the ultrasonic wave in F₂ and F₁₀ follows the same law, that is, as the bedding thickness increases, the propagation distance of the ultrasonic signal in F₁₀ is shorter in the same time interval, indicating that the ultrasonic wave velocity is smaller in F₁₀. In Fig. 9, different acoustic field amplitudes are used to characterize the ultrasonic transmission energy. It is obvious that the color of E₁₀ is lighter than E₂, indicating that the transmitted ultrasonic energy is smaller in E₁₀ as the vertical bedding thickness increases (Fig. 9b). Therefore, its attenuation coefficient is also larger (Fig. 8b). In the models with horizontal bedding F₂ and F₁₀, the same variation rule as E₂ and E₁₀ is revealed, that is, as the thickness of the bedding increases, the transmitted ultrasonic energy decreases (Fig. 9c, 9d), and the attenuation coefficient increases (Fig. 8b).

2.1.2. Influence of bedding density on ultrasonic characteristics

Under given bedding angle and thickness, the variation of bedding density was simulated by adjusting the number of bedding, that is, the number of bedding per unit area of the core section. With the vertical bedding thickness set at 0.4 mm, and the numbers of beddings were set at 0, 1, 2, 4, 5, 6, 10, 12, 15, 18, and 20, corresponding to the bedding density of 0, 0.08×10⁴, 0.16×10⁴, 0.32×10⁴, 0.40×10⁴, 0.48×10⁴, 0.80×10⁴, 0.96×10⁴, 1.20×10⁴, 1.44×10⁴, and 1.60×10⁴ beddings/m², respectively. The models of this group were coded G₁ to G₁₁. Models of another group with horizontal beddings were coded H1 to H₁₁, which have the bedding density same as the corresponding ones in the above group. The G₃, G₇ and H₃ and H₇ models which can better reflect the bedding density were selected to analyze the influence of the bedding density on the ultrasonic characteristics of the two groups of models, and analyze how the bedding density influences on the ultrasonic characteristics microscopically by wavefield snapshots. The bedding density of G₃ and G₇ are 0.16×10⁴ beddings/m² and 0.80×10⁴ beddings/m² respectively, that of H₃ and H₇ are 0.16×10⁴ beddings/m² and 0.80×10⁴ beddings/m² too.

The calculation results of wave velocity show that with the increase of bedding density, the wave velocity decreases linearly both in the horizontal bedding and in vertical bedding, and that in the vertical bedding decreases less (Fig. 10a). The numerical calculation results of the attenuation coefficient show that the attenuation coefficient increases linearly with the increase of the bedding density (Fig. 10b). The calculated results of dominant ultrasonic frequency and amplitude show that the dominant ultrasonic frequency and amplitude in the vertical bedding change irregularly with the increase of the bedding density, but in the horizontal bedding, with the increase of the bedding density and when the density is less than 0.3×10⁴ beddings/m², the dominant frequency and amplitude decrease, and when the density is more than 0.3×10⁴ beddings/m², the dominant frequency and amplitude rise slowly till basically keep stable (Fig. 10c).

The bedding density was preset in the models. Because the G₇ bedding density is larger than the G₃ bedding density, the H₇ bedding density is larger than the H3 bedding density, the ultrasonic wave propagated from the bottom to the top of the core, and the ultrasonic wave passed through more interfaces between the bedding and the skeleton in G₇ than in G₃. The snapshots of the wave field show that the ultrasonic signal propagates a shorter distance in G₇ than in G₃ in the same time interval (Fig. 11), so the ultrasonic velocity is also smaller in G₇ (Fig. 10a). The propagation of ultrasonic wave in horizontal bedding in H₇ and H₃ follows the same rule, that is, the propagation distance of the ultrasonic signal is shorter in H₇ than in H₃ in the same time interval, indicating that the ultrasonic wave velocity is also lower in H₇ than in H₃. In terms of attenuation, as the bedding density increases, G₇ has more beddings than G₃, and thus stronger heterogeneity than G₃. Due to the unequal modulus between the bedding and skeleton, a certain strain phase difference was generated between the skeleton and the bedding, so internal friction was generated at the interface, and the mechanical energy was converted into thermal energy, which greatly consumed energy^[19]. In Fig. 11, different acoustic field amplitudes are used to characterize the ultrasonic transmission energy. The G₇ color is lighter than G₃ (Fig. 11), therefore the ultrasonic transmission energy in G₇ is apparently lower than that in G₃, and the attenuation factor of G₇ is larger than that of G₃ (Fig. 10b). Similarly, H3 and H₇ show the same law, that is, as the bedding density increases, the ultrasonic transmission energy in H₇ is smaller, indicating that its attenuation coefficient is larger.

2.1.3. Influence of bedding angle on ultrasonic characteristics

With the bedding thickness and density kept constant, the models were set at the bedding angle of 0°, 10°, 20°, 30°, 40°, 50°, 60°, 70°, 80°, 90° and coded as I₁ to I₁₀ respectively. In this paper, I₃ and I₈ with obvious bedding angle features were selected to analyze the influence of bedding thickness on ultrasonic characteristics microscopically with their wave field snapshots. The bedding angle of I₃ and I₈ are 30° and 80° respectively.

The numerical simulation results show that as the bedding angle increases, the wave velocity decreases in a power function (Fig. 12a), while the attenuation coefficient increases linearly (Fig. 12b).

The bedding angle of the model I₈ preset was larger than the bedding angle of I₃, the ultrasonic wave propagated from the bottom to the top of the core, and under the premise of the same bedding density, I₈ has more bedding interfaces than I₃. Their wave field snapshots show that the distance that the ultrasonic signal can propagate in I₈ is shorter than that in I₃ in the same time interval (Fig. 13a, 13b), and the wave velocity in I₈ is lower (Fig. 12a). In terms of attenuation, since there are more interfaces between the skeleton and the bedding in I₈ than in I₃, there are more paths that hinder ultrasonic propagation in I₈ than in I₃. In Fig. 13, different acoustic field amplitudes are used to characterize the ultrasonic transmission energy. The color of I₈ is generally lighter than that of I₃, indicating that the transmitted ultrasonic energy of I₈ is lower, so the attenuation coefficient in I₈ is larger than in I₃ (Fig. 12b).

Based on the numerical simulation results of shale, the sensitivity of acoustic parameters to the bedding structure was quantitatively analyzed by the grey system correlation method^[20]. Firstly, acoustic parameters such as ultrasonic velocity and attenuation coefficient were used as reference sequences, and the structural parameters such as bedding angle, density and thickness were taken as comparison sequences:

The data was averaged dimensionless, i.e.

The correlation coefficient of the sequence is:

The correlation coefficient of the entire series is:

In the acoustic parameters, apart from the velocity that is sensitive to the bedding angle and density (with correlation value of over 0.80), the attenuation coefficient is also sensitive to the bedding thickness (with correlation value of 0.89), while the dominant amplitude and frequency are less correlated with the bedding structure (with correlation value of below 0.60) (Table 1). Therefore, in the process of utilizing the acoustic parameters to characterize the shale bedding structure in the future, in addition to the conventional sensitive factor, velocity, the attenuation coefficient can be employed jointly.

Using rock wave velocity to establish a dynamic mechanical parameter model is the method most commonly used in field development. The isotropic model and the lateral isotropic model are expressed as follows:

The isotropic mechanical parameter prediction model is:

The lateral isotropic mechanical parameter prediction model is^[21]:

It is found through numerical simulation that it’s not feasible to establish the dynamic mechanical parameter model of shale only with wave velocity, and the attenuation coefficient sensitive to the anisotropic characteristics should be incorporated to achieve better results. The wave velocity and attenuation coefficient of the Longmaxi Formation shale in the Pengshui area of southeastern Chongqing measured by experiments, and the McQuaid and global optimization algorithms^[22] were used to establish the dynamic mechanical parameter model considering dimension of the shale:

The density of the cores from Wells ZY1 YY1 were measured by physical properties experiments, the elastic modulus and Poisson’s ratio were obtained by mechanics experiments, and the P and S wave velocities were obtained by ultrasonic transmission experiments (Table 2). Then the models proposed in this paper were verified. The experimental results in Table 2 were taken into the isotropic model established by equations (14) and (15), and into the lateral isotropic model established by equations (16) and (17) to verify the models.

The results show that the dynamic and static mechanical parameters calculated by the models are in good agreement with experiment results. In Well ZY1, the correlation coefficient between static and dynamic Poisson’s ratios is up to 0.83, and the correlation coefficient between static and dynamic elastic moduli is 0.84. In Well YY1, the correlation coefficient between static and dynamic Poisson’s ratios is 0.82, and the correlation coefficient between static and dynamic elastic moduli is 0.81 (Fig. 14). The simulation effect of the new model proposed in this paper is better than the commonly used isotropic model and lateral isotropic model. Here the isotropic model refers to the medium with the same physical properties in different directions, and is a homogeneous model; the lateral isotropic model refers to the model with the same elasticity in any direction perpendicular to the axis, and is a heterogeneous model. The isotropic model only targets isotropic rocks, and does not consider the bedding development of shale. The lateral isotropic model regards shale as a laterally isotropic medium with structure too ideal. The model proposed in this paper is based on wave velocity and attenuation, the most sensitive parameters to shale bedding structure, and established by the general global optimization algorithm. The models can better reflect the influence of shale bedding structure, the dynamic moduli calculated by the models are generally larger than the corresponding static moduli, and the dynamic Poisson’s ratio to the static Poisson’s ratio is basically 1 to 1. They are basically consistent with the relationship obtained by Kuhlman et al.^[23].

The logging data of Well ZY2 in southeastern Chongqing was collected to carry out the inversion of the dynamic parameters of the rock mechanics, and the results of the static mechanical parameters were obtained by core experiments. Fig. 15 shows that the dynamic elastic modulus is generally larger than the static elastic modulus, and the dynamic Poisson’s ratio is very close to the static Poisson’s ratio, which is basically 1 to 1, in the inverted depth range. These are consistent with the model results mentioned above. Meanwhile, rock mechanics experiments were carried out on 10 cores from different depths of the Longmaxi Formation in Well ZY2, and inversion data of field logging data was verified by the experimental results.

Table 3 shows that the inverted mechanical parameters at 10 different depths are basically consistent with the experimental results of cores. Through calculation, the relative error of the Poisson’s ratio is less than 10%, and 6.94% on average; the relative error of the elastic modulus is 21.25% at maximum and 15.30% on average. Clearly, the mechanical profile inverted by the logging data can accurately reflect the real rock mechanical parameters. The dynamic mechanical parameter models established in this paper are suitable for the shale with developed bedding structure.

Based on the wave theory, a numerical simulation method of ultrasonic transmission experiment on shale was worked out by using the staggered grid finite difference method. The correlation coefficient between the simulated waveform and the experimental waveform is greater than 80%, which indicates that the numerical simulation method can effectively simulate the ultrasonic transmission experiment.

In the ultrasonic parameters, the acoustic velocity is a conventional sensitivity factor used to characterize shale bedding structure, but the normalized attenuation coefficient is needed to comprehensively describe the shale bedding to make the result more accurate.

Based on the wave velocity and attenuation coefficient of ultrasonic transmission experiments, the mechanical parameter models of the Lower Silurian Longmaxi Formation in southeastern Chongqing were established. The correlation between the dynamic and static parameters calculated by these models is better than the traditional model. The predicted values of the rock mechanics profile obtained from inverting logging data and the model are in good agreement with the experimental results.

A, B—discrete values of velocities u and ν, dimensionless;

D, E, F—discrete values of stresses τ_xx, τ_xz, τ_zz, dimensionless;

E₁, E₂—elastic moduli parallel and perpendicular to the isotropic plane, MPa;

i, j, k—time and space grids, dimensionless;

G₂—shear modulus in the isotropic plane, MPa;

p, q—parameters, dimensionless;

R₁—Acoustic displacement, m;

U, V—x and z velocities, m/s;

τ_xx, τ_xz, τ_zz—stress components, MPa;

Y_i, X_i—reference sequence and comparison sequence, dimensionless;

ξ_i(k)—correlation coefficient of X_i(k) and Y_i(k) at point k, dimensionless;

r_i—correlation degree of X_i(k) and Y_i(k), dimensionless;

ν₁, ν₂—Poisson’s ratios parallel and perpendicular to the isotropic plane, dimensionless;

V_p, V_s—P and S wave velocities, m/s;

λ, μ—Lame constant, dimensionless;

ρ_j—medium density, g/cm³;

ρ—core density, g/cm³;

χ—normalized attenuation coefficient, dimensionless.