Azimuthal elastic impedance-based Fourier coefficient variation with angle inversion for fracture weakness

Quantitative inversion of fracture weakness plays an important role in fracture prediction. Considering reservoirs with a set of vertical fractures as horizontal transversely isotropic media, the logarithmic normalized azimuthal elastic impedance (EI) is rewritten in terms of Fourier coefficients (FCs), the 90° ambiguity in the azimuth estimation of the symmetry axis is resolved by judging the sign of the second FC, and we choose the FCs with the highest sensitivity to fracture weakness and present a feasible inversion workflow for fracture weakness, which involves: (1) the inversion for azimuthal EI datasets from observed azimuthal angle gathers; (2) the prediction for the second FCs and azimuth of the symmetry axis from the estimated azimuthal EI datasets; and (3) the estimation of fracture weakness combining the extracted second FCs and azimuth of the symmetry axis iteratively, which is constrained utilizing the Cauchy sparse regularization and the low-frequency regularization in a Bayesian framework. Tests on synthetic and field data demonstrate that the 90° ambiguity in the azimuth estimation of the symmetry axis has been removed, and reliable fracture weakness can be obtained when the estimated azimuth of the symmetry axis deviates less than 30°, which can guide the prediction of fractured reservoirs.


Introduction
Subsurface fractures contribute to providing pathways for fluid flow and increasing the permeability of reservoirs. The key parameters of great interest for fractured reservoir exploration are the distribution and orientation of fracture systems, which play an important role in optimizing production from naturally fractured reservoirs (Sayers 2009;Far et al. 2014;Den Boer and Sayers 2018).
The seismic fracture prediction can be classified into two basic types: One is the seismic fracture qualitative prediction, which is mainly based on the discontinuity of fractures and uses the prestack seismic diffraction wave imaging approach or poststack geometric seismic attributes (such as curvature, coherence, and discontinuity) to qualitatively describe macroscale fractures, such as faults. The other is the seismic fracture quantitative prediction, which is mainly based on seismic azimuthal anisotropy (such as reflection amplitude, velocity, and attenuation) to predict mesoscale fractures, including obtaining fracture density and direction by ellipse fitting using azimuthal anisotropy attributes and quantitatively inverting fracture parameters using azimuthal prestack seismic data (Liu et al. 2015;Chen et al. 2016;Liu et al. 2018;Yuan et al. 2019).
Fracture detection utilizing seismic anisotropy has been an active field in the past decades. Hydrocarbon reservoirs with vertically aligned fractures can be described by horizontally transverse isotropic (HTI) media caused by a single set of rotationally invariant fractures in isotropic background Edited by Jie Hao 1 3 rock. The amplitude and velocity of seismic waves propagating in HTI media usually exhibit significant azimuthal anisotropy. Reflection amplitude is superior to seismic velocities in characterizing fractured reservoirs due to its higher vertical resolution and sensitivity to the properties of a reservoir (Far et al. 2013). Much work has been done in deriving linear P-wave reflection coefficients in HTI media (Rüger 1998;Pšencík and Martins 2001;Shaw and Sen 2006;Chen et al. 2014a;Pan et al. 2018a) and directly characterizing fractured zones using azimuthal P-wave reflectivity data (Mallick et al. 1998;Gray and Todorovic-Marinic 2004;Bachrach et al. 2009;Mahmoudian et al. 2015;Wang et al. 2019).
In recent years, amplitude variation with angle and azimuth (AVAZ) inversion combining with the linear slip theory (Schoenberg and Douma 1988;Schoenberg and Sayers 1995) used for modeling fractures embedded in isotropic host rock is widely applied to predict fractured reservoirs. Based on the rock physics model, Chen et al. (2014a) proposed the AVAZ inversion method to estimate elastic and fracture weakness parameters. Pan et al. (2019) presented Bayesian AVAZ direct inversion for fluid indicator and fracture weakness in an oil-bearing fractured reservoir. AVAZ inversion is an ill condition; a large number of unknown parameters and the cross talk between the elastic and fracture parameters inevitably make several-term AVAZ inversion unstable (Downton et al. 2006;Bachrach 2015). One tentative approach to reduce the number of unknown parameters is to implement azimuthal seismic amplitude difference inversion for fracture weakness estimation (Chen et al. 2017a;Pan et al. 2017;Xue et al. 2017). In addition, Downton andRoure (2011, 2015) rewrote azimuthal P-wave reflectivity in terms of Fourier coefficients (FCs) and proposed an azimuthal Fourier coefficient elastic inversion for fracture parameters estimation. Barone and Sen (2018) implemented the application of real seismic data targeting the Haynesville Shale using a Fourier azimuthal amplitude variation fracture characterization method. Nevertheless, the suitability and reliability of AVAZ inversion in the presence of noise are still controversial. Moreover, it is a challenging task to extract reasonable space variant wavelets of each incident angle at different azimuths used for AVAZ inversion. This can be addressed by azimuthal elastic impedance (EI) extended from the concept of EI (Connolly 1999;Whitcombe 2002;Wang et al. 2006;Yin et al. 2013;Zong et al. 2013Zong et al. , 2016Mozayan et al. 2018). In the last decades, there were abundant studies on the application of azimuthal EI in anisotropic inversion (Martins 2006;Chen et al. 2014bChen et al. , 2017bPan et al. 2018b;Pan and Zhang 2019).
Our study is the extension of Downton and Roure (2015) research, which stabilizes the inversion process by treating the fracture weakness separately from the elastic parameters using FCs. Different from traditional fracture weakness inversion using azimuthal reflection amplitude or azimuthal EI or azimuthal reflection amplitude FCs, in this paper, a novel azimuthal EI-based FC variation with angle inversion method is proposed. The normalized azimuthal EI equation in HTI media is rewritten by using the Fourier series expansion method, the azimuth of the symmetry axis is estimated without 90° ambiguity using FCs, and FCs with the highest sensitivity to fracture weakness are chosen to estimate fracture weakness, which is constrained utilizing the Cauchy sparse regularization and the low-frequency regularization in a Bayesian framework. The influence of the azimuth of the symmetry axis on fracture weakness inversion is analyzed, and the proposed approach is demonstrated on both synthetic and real data.

Fourier coefficients of normalized azimuthal EI in logarithm domain
The P-wave reflection coefficient in HTI media is given by (Pan et al. 2017): with where , , and represent P-wave velocity, S-wave velocity, and density of the background isotropic medium, respectively. The superscriptrepresents the average of elastic parameters. g = 2 2 is the ratio of the squared vertical S-and P-wave velocities in the background isotropic medium, and Δ N and Δ T are the normal and tangential fracture weakness of layers. If the HTI model is induced by penny-shaped cracks, Δ T gives a direct estimate of the crack density and the ratio Δ N ∕ Δ T is a sensitive indicator of fluid saturation (Bakulin et al. 2000). Δ Δ N represents the difference in normal weakness between the upper and lower layers, and Δ Δ T represents the difference in tangential weakness between the upper and lower layers. is the incident angle, D( , ) = −2g cos 2 sin 2 + sin 2 cos 2 sin 2 tan 2 (1 − 2g) + cos 4 sin 2 tan 2 (1 − g) , E( , ) = 2g cos 2 sin 2 1 − sin 2 tan 2 , 1 3 and = − sym is the angle between the observed azimuth and azimuth sym of the symmetry axis. Pan et al. (2017) further derived the normalized azimuthal EI equation as follows: where the subscript zero represents the constant medium properties of elastic parameters. The dimensionality dependence on the incident angle is addressed by the introduction of EI 0 = (Whitcombe 2002). Taking logarithms on both sides of Eq. (2) yields: with Equation (3) is the basis of conventional EI variation with angle and azimuth inversion for elastic and fracture weakness parameters. However, the robustness of inverting fracture weakness from azimuthal EI is still controversial due to the coupling between elastic and fracture weakness parameters. Downton andRoure (2011, 2015) rewrote the azimuthal P-wave reflection coefficient using FCs. We can also describe the normalized azimuthal EI in logarithm domain in terms of FCs and rearrange Eq. (3) as follows: with where r 0 ( ) represents the DC component FC, which involves the three-term AVO expression with the addition of fracture weakness terms, r 2 ( ) and r 4 ( ) represent the second (4) L EI ( , ) = r 0 ( ) + r 2 ( ) cos 2 − sym + r 4 ( ) cos 4 − sym (5) g − 1 sin 2 tan 2 , G( ) = g sin 2 − 1 4 g sin 2 tan 2 , H( ) = g(2g − 1) sin 2 + g(g − 1) sin 2 tan 2 , I( ) = g sin 2 , J( ) = − 1 4 g 2 sin 2 tan 2 , K( ) = − 1 4 g sin 2 tan 2 , and fourth FCs, respectively, which are only related to fracture weaknesses and independent of elastic parameters, and the phase of the second and fourth FCs is related to azimuth of the symmetry axis. Equation (4) can be further written as the weighted sum of cosine and sine waves by: where a n ( ) and b n ( ) (n = 0, 2, 4) represent the nth FCs, respectively. For the case of X regularly sampled data, a discrete Fourier transform (DFT) can be used to calculate a n ( ) and b n ( ): Note that Nyquist criterion must be met so that a n ( ) and b n ( ) can be obtained by DFT. In order to estimate a 2 ( ) and b 2 ( ) , the data must be sampled no coarser than 45°. In order to estimate a 4 ( ) and b 4 ( ) , the data must be sampled no coarser than 22.5°.

Sensitivity analysis of the FCs
The second and fourth FCs are related to fracture weaknesses and independent of elastic parameters; we use the isotropic medium model parameters given by Rüger and Tsvankin (1997) and change the normal and tangential weaknesses, respectively, to analyze the sensitivity of the second and fourth FCs to normal and tangential weaknesses. The P-wave velocity, S-wave velocity, and density of the isotropic medium model are 4.500 km/s, 2.530 km/s, and 2.800 g/cm 3 , respectively. The normal and tangential weaknesses increase from 0 to 0.2, respectively, in increments of 0.05. We use Eqs. (6) and (7) to calculate the second and (8) fourth FCs. Figure 1a, b shows the second FC r 2 ( ) variation with the normal and tangential fracture weaknesses. Figure 2a, b shows the fourth FC r 4 ( ) variation with the normal and tangential fracture weaknesses. It can be seen that when the incident angle is less than 30°, r 4 ( ) hardly changes with normal and tangential weaknesses, and r 2 ( ) changes significantly with normal and tangential weaknesses. Compared with r 4 ( ) , r 2 ( ) is more sensitive to normal and tangential weaknesses.
Since a 2 ( ) = r 2 ( ) cos 2 sym and b 2 ( ) = r 2 ( ) sin 2 sym , a 2 ( ) is a scaled version of r 2 ( ) by the scalar cos 2 sym and b 2 ( ) is a scaled version of r 2 ( ) by the scalar sin 2 sym . In the case when one of a 2 ( ) and b 2 ( ) is zero, the other one is r 2 ( ) or −r 2 ( ) . As shown in Figs. 1a, b and 3a, b, both r 2 ( ) and −r 2 ( ) are sensitive to fracture weaknesses. As a result, there is always one in a 2 ( ) and b 2 ( ) which is sensitive to fracture weaknesses. In the case when a 2 ( ) and b 2 ( ) are the same, both are to fracture weaknesses, so both a 2 ( ) and b 2 ( ) are sensitive to fracture weaknesses. In other cases, at least one of a 2 ( ) and b 2 ( ) is sensitive to fracture weaknesses, and its sensitivity to fracture weaknesses is greater than the sensitivity of √ 2r 2 ( ) quently, no matter how sym changes, at least one of a 2 ( ) and b 2 ( ) is sensitive to fracture weaknesses, which indicates that a 2 ( ) and b 2 ( ) can be simultaneously combined to invert for fracture weaknesses.

Resolving the 90° ambiguity in the azimuth estimation of the symmetry axis
According to Eqs. (10) and (11), a 2 ( ) and b 2 ( ) are functions of sym . In order to estimate fracture weaknesses from a 2 ( ) and b 2 ( ) , sym should be estimated first, which can be obtained by: Since the signs of r 2 ( ) and r 4 ( ) cannot be determined, there is a 90° ambiguity in sym estimate using a 2 ( ) and b 2 ( ) , and there is a 45° ambiguity in sym estimate using a 4 ( ) and b 4 ( ) ; we propose an approach to resolve the ambiguity. Firstly, sym is directly calculated from a 2 ( ) and b 2 ( ) by: Since a 2 ( ) = r 2 ( ) cos 2 sym , b 2 ( ) = r 2 ( ) sin 2 sym , for the case of r 2 ( ) > 0 , if the sign of a 2 ( ) is the same as the sign of cos 2 sym , and the sign of b 2 ( ) is the same as the sign of sin 2 sym , the estimated sym is correct, otherwise the estimated sym = sym + 90 • . For the case of r 2 ( ) < 0 , if the sign of a 2 ( ) is opposite to the sign of cos 2 sym , and the sign of b 2 ( ) is opposite to the sign of sin 2 sym , the estimated sym is correct, otherwise the estimated sym = sym + 90 • .
Combining Eqs. (6), (21), and (22) gives: Thus, for dry (gas-filled) cracks, r 2 ( ) is related to the incident angle, crack density, and the ratio g , and since the crack density e is always positive, the sign of r 2 ( ) is only related to and g . As shown in Fig. 4a-d, for the case of different incident angles, r 2 ( ) increases with the ratio g and changes sign with increasing g , and the slope of r 2 ( ) gradually increases with crack density. However, the change of the crack density does not change the corresponding ratio g for r 2 ( ) changing sign, while the change of the incident angle changes the corresponding ratio g for r 2 ( ) changing sign, and the corresponding ratio g for r 2 ( ) changing sign becomes larger with increasing incident angle. Note that for the case of incident angle ≤ 30 (in Fig. 4a- To judge the sign of r 2 ( ) using the above method requires knowing fracture infill (content) in advance, but fracture infill is usually unknown, so anisotropic rock physical model Chen et al. 2014a;Pan et al. 2017) can be used to estimate fracture weaknesses, then calculate r 2 ( ) , and judge its sign to further resolve ambiguity.

Azimuthal EI-based FC variation with angle inversion for fracture weakness
For each incident angle, the normalized azimuthal EI in logarithm domain can be decomposed into a series of FCs, the second and fourth FCs are only related to fracture weakness and independent of elastic parameters, and the previous analysis suggests that we use second FCs to estimate fracture weakness.
In the case of M incident angle and N time sample, combining Eqs. (10) and (11), we can get the following matrix expression: .
Here, the superscript T denotes the transpose of a matrix, and the symbol diag denotes a diagonal matrix; the subscripts m represents the mth angle of incidence.
Equation (24) can be expressed in matrix form as: Based on Bayesian theory (Buland and More 2003), the posterior probability density function of the model parameters can be expressed as: We assume that the likelihood function p( | ) obeys a Gaussian probability distribution, which is given by: where 2 n is the variance of noise. Cauchy probability distribution can suppress random noise and produce a sparse solution with high vertical resolution (Downton 2005;Chen et al. 2017b;Yuan et al. 2018); we assume that the prior probability distribution p( ) obeys Cauchy probability distribution as follows: where 2 m is the variance of model parameters. Combining Eqs. (26), (27), and (28), we can get the initial objective function for the MAP solution inversion, which is expressed as: We combine the low-frequency information constraint obtained from well logs; the final objective function is obtained as follows: where Δ N and Δ T are weighting factors for the normal and tangential fracture weaknesses, respectively. and represent initial models for the normal and tangential fracture weaknesses, which can be obtained by using the estimation results from the anisotropic rock physics model Chen et al. 2014a;Pan et al. 2017). The iteratively reweighted least squares (IRLS) method is employed to solve Eq. (30) for estimating the normal and tangential fracture weaknesses.

Workflow of azimuthal EI-based FC variation with angle inversion for fracture weakness
A workflow of azimuthal EI-based FC variation with angle inversion is proposed to estimate fracture weaknesses; we conclude the whole process as follows: 1. The constrained sparse spike inversion (CSSI) for azimuthal EI datasets from azimuthal angle stacks. The inputs include azimuthal angle stacks, azimuthal seismic wavelets, and low-frequency azimuthal EI models of each incident angle. Since azimuth sym of the symmetry axis is unknown, low-frequency azimuthal EI models are obtained by interpolating isotropic EI calculated from well log data. 2. The prediction for a 2 ( ) and b 2 ( ) from all estimated azimuthal EI datasets and estimation for sym from a 2 ( ) and b 2 ( ) . The ambiguity in sym estimate is resolved by judging the sign of r 2 ( ) , which is obtained by combining fracture weaknesses from the anisotropic rock physics model Chen et al. 2014a;Pan et al. 2017) and elastic parameters from well log.
Note that for the case of regularly sampled azimuthal EI datasets, a 2 ( ) and b 2 ( ) can be calculated using a DFT shown in Eqs. (14) and (15). For the more complex case of irregularly sampled EI datasets in azimuth, a least squares inversion may be performed by using Eq. (8). 3. The estimation for fracture weaknesses using the extracted a 2 ( ) , b 2 ( ) , and sym iteratively, which is constrained utilizing the Cauchy sparse regularization and the low-frequency regularization in a Bayesian framework (Pan and Zhang 2019).

Synthetic examples
In order to validate the proposed approach, a five-layer medium model is established (as indicated in Table 1), in which the first, third, and fifth layers are isotropic medium, and the second and fourth layers are HTI medium resulting from dry cracks; the model parameters are given by Rüger and Tsvankin (1997). Note that the fracture symmetry axis azimuths of the second and fourth layers are assumed to be 60° and 120°, respectively. Thomsen anisotropy parameters are given in Table 1; the normal and tangential weaknesses are obtained using the approximation given by Bakulin et al. (2000): Synthetic data are generated from model parameters given in Table 1 convolved with a 35 Hz Ricker wavelet, including four azimuthal angle gathers (the given azimuths are 0°, 45°, 90°, and 135°, respectively), and the incident angle ranges from 5° to 35°. Firstly, angle gathers for each azimuth are stacked to obtain azimuthal angle stacks, including three average incident angles (small incident angle 10°, stacked by 5°-15°; middle incident angle 20°, stacked by 15°-25°; large incident angle 30°, stacked by 25°-35°); for each azimuth, the azimuthal EI datasets are obtained from azimuthal angle stacks and used to calculate a 2 ( ) and b 2 ( ) . Since a 2 ( ) and b 2 ( ) are related to sym and the ratio g , the ratio g can be calculated from a constant or slowly varying prior known background model (Buland and More 2003). In order to analyze the effect of sym on fracture weakness inversion, we combine a 2 ( ) , b 2 ( ) , and different sym (the azimuth deviation Δ sym of the symmetry axis is 0°, 15°, 30°, 45°, 60°, 75°, and 90°, respectively) to estimate fracture weaknesses. Figure 5a shows the synthetic azimuthal angle gathers without noise; Fig. 6 shows the inversion results of fracture weaknesses in different Δ sym cases, the black curve represents the true models, and the dashed lines of different colors, respectively, represent the inversion results in the case of different Δ sym . It can be seen that when Δ sym is 0°, the inversion results are completely consistent with the true models. When Δ sym is 15°, the inversion results are slightly deviated from the true models. When Δ sym is 30°, the inversion results are more deviated from the true models, but they are basically consistent with the true model trends. As Δ sym continues to increase ( Δ sym is 45°, 60°, 75°, and 90°, respectively), the inversion results are gradually opposite to the true model trends and are no longer reliable, which indicates that fracture weakness inversion results obtained by the proposed approach is reliable when Δ sym is less than 30°. sym should be estimated before fracture weakness inversion from a 2 ( ) and b 2 ( ) ; we propose a stable approach to estimate sym without ambiguity and also use the above model to validate the proposed approach. Figure 7a shows the directly estimated sym , and Fig. 7b shows the corrected sym using the proposed approach. Note that the directly estimated sym of the second layer (pink rectangle) and fourth layer (green rectangle) is 150° and 30°, respectively, both having a 90° ambiguity, while the relevant corrected sym is 60° and 120°, respectively, which indicates that the 90° ambiguity is successfully resolved using the proposed approach. It is noted that in the proposed approach, the sign of r 2 ( ) should be determined first. Here, r 2 ( ) is calculated by smoothing model parameters for 100 times. In order to show the performance of the proposed approach in synthetic noisy data, white Gaussian noise with different signal-to-noise (S/N) ratios is added to the resulting synthetics; the relevant synthetic noisy results are shown in Fig. 5b, c, respectively. The relevant directly estimated sym is shown in Figs. 8a and 9a, respectively, and the relevant corrected sym is shown in Figs. 8b and 9b, respectively. As can be seen from Fig. 8a, b, for the case of S/N = 5, in the second layer, the sym directly estimated at small, middle, and large incident angles is approximately 135°, 143°, and 149°, respectively, and the relevant corrected sym is approximately 45°, 53°, and 59°, respectively. In the fourth layer, the sym directly estimated at small, middle, and large incident angles is approximately 173°, 27°, and 28°, respectively, and the relevant corrected sym is approximately 173°, 117°, and 118°, respectively. As can be seen from Fig. 9a, b, for the case of S/N = 2, in the second layer, the sym directly estimated at small, middle, and large incident angles is approximately 25°, 141°, and 151°, respectively, and the relevant corrected sym is approximately 25°, 51°, and 61°, respectively. In the fourth layer, sym directly estimated at small, middle, and large incident angles is approximately 171°, 23°, and 35°, respectively, and the relevant corrected sym is approximately 171°, 113°, and 125°, respectively. It can be seen that sym directly estimated at small incident angles is unreliable, and sym directly estimated at middle and large incident angles is more stable and reliable but has a 90° ambiguity, while the 90° ambiguity is successfully resolved using the proposed approach. It indicates that sym can be estimated by using seismic data at middle and large incident angles. If there are multiple sym estimated at middle and large incident angles, the average of multiple sym can be a more stable and reliable result.
The fracture weakness can be obtained by combining a 2 ( ) , b 2 ( ) , and the estimated sym . The above model is also used to validate the proposed fracture weakness inversion approach. Here, the average of sym estimated at middle and large incident angles is taken as the final sym for fracture weakness inversion. Figure 10a-c shows the comparisons between the inversion results (red) and the true models (black) of the normal and tangential weaknesses. Here, the initial models (green) of fracture weaknesses are Table 1 Elastic and anisotropic parameters of the five-layer model (Rüger and Tsvankin 1997) Layers , km/s , km/s , g/cm 3 the smooth results of true models. It can be seen that the inversion results are basically consistent with the true models in the case of no noise. In the case of S/N being 5 and 2, the inversion results still show good consistency with the true models, which further verifies the robustness of the proposed approach.

Field data example
A 2D line through single well acquired in western China is used to further verify the stability of the proposed approach. Noise and multiples must be carefully removed to get the true-amplitude-processed gathers before inversion. When sorting azimuths, we should ensure sufficient coverage and S/N at each azimuth division. Here, the data were sorted into four azimuth sectors, including 0° (stacked by − 22.5°-22.5°), 45° (stacked by 22.5°-67.5°), 90° (stacked by 67.5°-112.5°), and 135° (stacked by 112.5°-157.5°). For each divided azimuth, seismic data is stacked over the incident angle to obtain three angle stacks, including small angle 10° (stacked by 5°-15°), middle angle 20° (stacked by 15°-25°), and large angle 30° (stacked by 25°-35°). We first implement CSSI for azimuthal EI datasets from azimuthal angle stacks. Figure 11a-c shows the azimuthal angle stacks at different incident angles. Figure 12a-c shows the corresponding azimuthal EI inversion results. We next implement DFT for a 2 ( ) and b 2 ( ) from all estimated azimuthal EI datasets and estimation for sym from a 2 ( ) and b 2 ( ) . Under the HTI medium assumption, the fast P-wave orientation is usually parallel to the fracture strike and perpendicular to the fracture symmetry axis (Delbecq et al. 2013). Since there is fast P-wave orientation data in the work area, we use the fast P-wave orientation to rotate 90° to calculate sym , which is used as the true sym for comparison with the sym estimated by the proposed approach. Figure 13a shows the true sym calculated from the fast P-wave orientation. Figure 13b-d shows the estimated sym at different incident angles. It can be seen that the estimated sym at small incident angle is the worst and the estimated sym at middle and large incident angles is relatively better. The estimated sym near the well location is better than that away from the well location, which may be caused due to the unreasonable fracture weakness models away from the well location. The average of estimated sym at middle and large incident angles is taken as the final estimated sym . Figure 14a shows the final estimated sym , and Fig. 14b shows the absolute values of the error between the final estimated sym and true sym . It can be seen that most of the final estimated sym is basically consistent with true values, and most of the absolute values of the error between the final estimated sym and true sym are less than 30°.
Finally, combining a 2 ( ) and b 2 ( ) , the fracture weakness inversion is performed using true sym and the final estimated sym , respectively. Figure 15a, b shows the inversion results of fracture weakness obtained using true sym ; Fig. 16a, b shows the inversion results of fracture weakness obtained using the final estimated sym . It can be seen that the inversion results near the well location are basically the same in both cases, and there are some differences in the inversion results away from the well location, which may be caused due to the error in the sym estimation. As a whole, the inversion results in both cases are basically the same. We plot the well logs for comparison with the inversion results, which show good agreement with the observation in wells. The large values of the inverted normal and tangential weaknesses near the well are consistent with the fractured zones in the reservoir. Therefore, the application of field data further verifies the validity and feasibility of the proposed approach.

Conclusions
Fractured reservoirs with a set of vertical fractures are equivalent to HTI medium under long-wavelength assumptions. A logarithmic normalized azimuthal EI equation is rewritten in terms of FCs; a novel azimuthal EI-based FC variation with angle inversion is proposed to estimate fracture weakness. Finally, the synthetic and real data are used to verify ig. 7 Comparisons between the directly estimated sym and the relevant corrected sym at different incident angles (circle lines in different colors) for the case of no noise. a The directly estimated sym , and b the relevant corrected sym . Note that the true sym of the second layer (pink rectangle) and the fourth layer (green rectangle) is 60° and 120°, respectively the feasibility and rationality of the proposed approach, and the following conclusions are obtained: 1. The approach based on reflection amplitude for estimating azimuth of the symmetry axis has a 90° ambiguity. In this paper, the ambiguity is eliminated by judging the sign of the second FC r 2 ( ) . If fracture infill is known, for example, for fluid-filled cracks, there is always r 2 ( ) > 0 . If the cracks are dry (gas-filled), in this case it becomes complicated to determine the sign of r 2 ( ) , which is related to incident angle and the ratio g . For the case of ≤ 30 • , if g > 0.4 , there is always r 2 ( ) > 0 . For the case of > 30 • , if g < 0.4 , there is always r 2 ( ) < 0 . While fracture infill is usually unknown, anisotropic rock physical model can be used to estimate fracture weakness, then calculate r 2 ( ) , and judge its sign to fur-ther resolve ambiguity. In addition, the analysis shows that seismic data at middle and large incident angles should be used to estimate sym . In the case that there are multiple sym estimated at middle and large incident angles, the average of multiple sym can be a more stable and reliable estimate, which contributes to reducing the uncertainty of the sym estimate. 2. The sensitivity analysis of the FCs to fracture weaknesses shows that the second FC is more sensitive to fracture weaknesses than the fourth FC; the second FCs obtained from azimuthal EI data can be simultaneously combined to estimate fracture weaknesses. However, azimuth of the symmetry axis should be estimated prior to the fracture weakness inversion, and the fracture weakness inversion is affected by the estimated azimuth of the symmetry axis. Tests on the synthetic and real 1 3 data show that the inversion results of fracture weaknesses obtained by the proposed approach are stable and reliable when azimuth of the symmetry axis deviates less than 30°. 3. The proposed approach for estimating azimuth of the symmetry axis and fracture weaknesses depends on the constructed anisotropic rock physical model. The correctness of the model affects the accuracy of the azimuth of the symmetry axis and fracture weakness estimation, so questions remain about the reliability of this inversion, and more work needs to be done to improve the reliability of fracture prediction.    Trace number EI, kg·m -2 ·s -1 ×10 6

Fig. 12
Inverted azimuthal EI data at small, middle, and large incident angles, in which a shows an average angle of 10° (small incident angle), b shows an average angle of 20° (middle incident angle), and c shows an average angle of 30° (large incident angle). The red curve represents the P-wave impedance at the well location Trace number ∆ N ∆ T Fig. 16 Inversion results of fracture weakness obtained using the final estimated sym , a the inverted normal weakness and b the inverted tangential weakness