Multiple elastic attribute fusion approach based on affinity propagation clustering strategy for gas hydrate reservoir identification

Pre-stack seismic inversion based on sensitive elastic parameters is critical in reservoir lithology prediction and geofluid identification. The ability of a single elastic attribute to identify a reservoir depends on its sufficient sensitivity to distinguish the target reservoir from the surrounding sediments. In general, high-dimensional data space composed of multiple elastic attributes is more conducive to describing reservoir characteristics. Therefore, a multiple elastic attribute fusion method using affinity propagation clustering strategy for gas hydrate reservoir identification is proposed. Rock-physics modeling is the most effective tool to determine the influence of microscopic physical parameters on macroscopic elastic response and to quantitatively evaluate the sensitivity of elastic attributes. Consequently, a rock-physics model of hydrate sediments considering the non-negligible shear properties of hydrates is constructed to clarify reservoir-sensitive elastic parameters. Additionally, a clustering evaluation indicator is defined to determine the optimal data clustering dimension in terms of feasibility and economy of the proposed method, and to avoid bias in the results due to data redundancy. It is shown that the 3D elastic attribute space consisting of shear modulus, Young's modulus, and S-wave velocity has the best discrimination ability for hydrate reservoirs. The logging data are used to verify the feasibility and effectiveness of the proposed method. Finally, the hydrate reservoir development is accurately discriminated by using the multiple elastic attributes yield from the pre-stack seismic inversion and combined with the fusion strategy.


Introduction
Gas hydrate is a new type of clean energy with huge reserves, which has become a research hotspot in exploration and development in recent years (Riedel et al. 2022). The reservoir identification method based on pre-stack seismic inversion aims to distinguish favorable targets from other sedimentary layers with the help of elastic attributes that have significant identification ability. According to the geological characteristics of different research areas and the properties of target reservoirs, a quantitative assessment of the sensitivity of the elastic parameters is first required. In addition, from the perspective of data dimensionality, the ability of a single attribute to describe the reservoir properties is always insufficient, even if this attribute is relatively more sensitive. Therefore, this study focuses on hydrate reservoirs, quantitatively evaluates the elastic parameters by constructing a rock-physics model, and uses a clustering strategy to fuse the select elastic attributes to achieve a comprehensive portrayal of the hydrate reservoir.
For hydrocarbon-bearing reservoirs, elastic parameters that more directly reflect lithology or geofluid are discovered and proposed, such as elastic impedance, Lamé parameter, Poisson's ratio, Young's modulus, bulk modulus, and Gassmann's fluid term, etc. (Goodway et al. 1997;Gray et al. 1999;Russell et al. 2003Russell et al. , 2011Wang et al. 2006;Zhang et al. 2009;Jones and Terry 2012;Zong et al. 2013Zong et al. , 2016Othman et al. 2021). However, there are still uncertainties and ambiguities in the quantitative identification of hydrate reservoirs affected by geological and geophysical factors and specific hydrate microscopic microstructures. The multiattribute clustering technique provides theoretical support for the integrated use of valuable information from seismic data for the accurate identification of target reservoirs (Singh et al. 2020). The K-means algorithm achieves data classification by determining the minimum Euclidean distance between each data point and the centroid of that cluster, and its simple structure is widely used in data mining (Mac-Queen 1967;Coléou et al. 2003;Giovanni and Lo 2018). It should be noted that the K-means and the corresponding improved algorithms require a given number of clusters using a priori information (Yang and Wang 2008). Principal component analysis (PCA) and self-organizing maps (SOM) are completely data-driven unsupervised clustering strategies with excellent applications in complex high-dimensional data analysis and seismic facies identification (Dumay and Fournier 1988;Matos et al. 2007;Duan et al. 2019). Fuzzy c-means (FCM) are also used for multiparametric analysis of geophysical methods for target detection, enhancing the robustness of clustering as it allows the same data to be divided into different subsets (Sun and Li 2015). Clustering strategy has been proven to be feasible for lithology discrimination in logging data with more satisfactory results Ojha 2021, 2022). However, the high correlation and nonlinearity of various seismic elastic attributes lead to less prominent category characteristics, which increases the difficulty of clustering implementation. Therefore, clustering methods with higher accuracy in elastic attribute clustering need to be explored. The affinity propagation (AP) clustering treats each data point equally, strictly complies with data distribution patterns, and efficiently identifies typical cluster members to give the best clustering schemes (Frey and Dueck 2007;Wang et al. 2013). Further, the AP strategy searches for a cluster center in the dataset through an information transfer mechanism, and that cluster center is the exact data point in the data to be clustered, which is more in line with the characteristics of spatial continuity of elastic parameters.
Rock-physics modeling is a means to link the microscopic pore structure, physical parameters, fluid properties, and macroscopic elastic properties of rocks (Mavko et al. 2009). Many rock-physics models suitable for hydrate layers have been proposed. Most attempts to better consider the pore structures, mineral types, pressure, and other geological conditions are based on effective medium theory. According to the hydrate distribution in sediments with different morphologies, the hydrate rock-physics model can be summarized as pore-filling model, load-bearing model and cementation model based on isotropic assumption, and fracture-filling model based on the anisotropic assumption (Helgerud et al. 1999;Ecker et al. 2000;Jenkins et al. 2005;Waite et al. 2009;Jana et al. 2015;Liu et al. 2018;Ghosh et al. 2021). Besides, considering the condition of unconsolidated rocks in hydrate-bearing sediments, Hertz-Mindlin theory is frequently referenced as the underlying theory in terms of rockphysics modeling (Mindlin 1949;Sain et al. 2010;Darrell and Camelia 2018;Pan et al. 2020). A reasonable hydrate rock-physics model clarifies the microscopic dominant factors affecting the elastic response of rocks, while quantitatively evaluating the sensitivity of elastic parameters to guide the selection of properties for clustering strategies.
In this study, we try to use a clustering strategy to integrate multiple elastic attributes, and the fused attributes can better reflect the reservoir characteristics. Firstly, a rockphysics model considering the shear properties of hydrates and the unconsolidated state of the rock is constructed to quantitatively assess the sensitivity of the elastic parameters of hydrate reservoirs. Secondly, the elastic attributes are selected as the initial data for clustering. Theoretically, the higher the dimensionality of the data, the more accurate the description of the reservoir. However, from the viewpoint of feasibility and economy, the more attributes, the workload increases exponentially. Therefore, the indicator of the clustering optimal dimension is defined by evaluating the intercluster distance concerning the intra-cluster data compactness. Finally, using the advantages of affinity propagation clustering strategy in efficiency and accuracy, the clustering results can clearly identify the hydrate reservoir.

Rock-physics modeling of hydrate reservoirs
A rock-physics model considering the solid-like properties of hydrates and the unconsolidated state of the rock by combining the applicability of the generalized Gassmann equation when the pores are filled with solids is constructed.

Solid matrix
The mineral composition of hydrate reservoirs is dominated by clastic minerals. In calculating the effective elastic modulus of the solid matrix, the Hashin-Shtrikman bound provides theoretically reasonable estimates. The general form of the Hashin-Shtrikman boundary in the case of multiple mineral mixtures is (Berryman 1995) where ± = ± (9K ± + 8 ± ) 6(K ± + 2 ± ) −1 ,K + and K − are the maximum and minimum bulk modulus of each solid components, + and − are the maximum and minimum shear modulus, and K i , i and f i are the bulk modulus, shear modulus, and volume content of the ith mineral component, respectively, N is the types contained in the mineral components. The elastic modulus of the solid matrix can be further obtained by the Hill average,

Dry frame
The addition of pores to the solid matrix forms the dry skeleton. Generally, we make assumptions about the pore shape and type, but the reasonableness of the assumptions is debatable due to the complexity of the hydrate microstructure, so the consolidation parameter model is used to calculate the elastic modulus of the dry frame (Lee 2008), (3) K ma = K HS + +K HS− 2 ma = HS + + HS− 2 . (4) where K d and d are bulk and shear modulus of dry frame, is porosity, and is consolidation parameter. The consolidation parameter model can describe the unconsolidated situation of hydrate-bearing sediments, and because it contains only one parameter to be estimated, it avoids the problem of difficult determination of parameters such as pore type and shape similar to those in DEM and other models.

Equivalent fluid
When the pore space of rock is filled with a non-viscous fluid such as water or gas, the shear modulus of the equivalent fluid is not considered. However, the shear modulus of hydrates is not negligible due to their solid-like properties. In addition, the investigation of the elastic properties of hydrates in different regions confirms the necessary consideration of the shear modulus of hydrates in the rock-physics modeling process (Table 1). The elastic modulus of equivalent fluid is calculated as, where K f and f are the bulk and shear modulus of equivalent fluid, S w is water saturation, K h and h are the bulk and shear modulus of gas hydrate, and K w is bulk modulus of water, respectively.

Saturated rock
For the solids filling the pore space of the sediment, Ciz and Shapiro (2007) proposed an elastic tensor equation for solidsaturated porous rocks, where S s , S d , S m and S if are the compliance tensor of saturated patchy, dry patchy, matrix and pore-filling material, respectively, S is the compliance tensor related to pore space of the dry frame. For the rock frame satisfies the isotropic assumption, the Eq. (9) can be expressed in the form of the elastic modulus, where: K sat and sat are the bulk and shear modulus of saturated rock, K dry and dry are the bulk and shear modulus of rock skeleton, K gr and gr are the bulk and shear modulus of matrix, K and are the bulk and shear modulus related to pore space, K if and if are the bulk and shear modulus related to filler inside pore space. Using the measured P-wave velocity in combination with the Monte Carlo (MC) algorithm, the consolidation parameters can be predicted to obtain more reasonable elastic calculated results, where E( ) is the error between the logging velocity V M Log and the model predicted velocity V C mod el , K , U , , are the bulk modulus, shear model, porosity and density of saturated rock during the model calculation, f gh , S gh are the volume content and saturation of hydrates, respectively. The solution of the objective function is the optimal estimate of consolidation parameter . (10)

Affinity propagation clustering strategy
AP clustering algorithm is a kind of clustering algorithm based on "information transfer" between data points. The basic principle is: firstly, the information transfer mechanism is used to search the clustering center in the data set and the affiliation between the data points and the data center; secondly, the data set to be clustered is divided into several subsets with specific meaning according to the affiliation. The advantages of AP algorithm include: the algorithm does not need to explicitly determine the parameters related to the number of clusters in the clustering process; the AP algorithm looks for "exemplars", i.e., clustering centroids that are actual points in the data set, as representatives of each class; the input to the algorithm can be either a symmetric similarity matrix or an asymmetric similarity matrix.
x 1 , x 2 , ⋯ , x n is the data sample set, and s is a matrix that inscribes the similarity between data points such that s(i, j) > s(i, k) when and only if the degree of similarity between x i and x j is greater than its similarity to x k . For the AP algorithm, all data points are potential clustering centers by alternating two message-passing to update two matrices: responsibility and availability matrices . The responsibility matrices describes the degree to which data object k fits as a clustering center for data object i , expressed as a message from i to k , defined as r(i, k) . Conversely, the availability matrices describes the degree of fitness of data object i to select data object k as its according clustering center, indicating the messages from k to i , defined as Step 0: Initialization definition of parameters. In the initial case, both the responsibility and availability matrices are defined as 0 matrices, Since all data points are potential clustering centers, the value s(i, j) on the diagonal of the similarity matrix s is assigned uniformly, which is generally the median or average of all similarities, called the preference. The algorithm can control the granularity of the clustering results by adjusting the preference, and the larger the    value of the preference, the greater the number of clustering results, and vice versa.
Step 1: Update responsibility matrix , Step 2: Update availability matrix , Step 3: To avoid oscillations, an attenuation factor is introduced when updating messages. Each message is set to be times the updated value of its previous iteration plus 1 − times the updated value of the current message, where ∈ (0, 1).
The above steps are iterated and the algorithm ends if these decisions remain unchanged after several iterations or if the algorithm executes more than the set number of iterations, or if the decisions about sample points in a small region remain unchanged after several iterations. The affinity propagation clustering algorithm is a data-driven unsupervised clustering strategy that does not require any prior information to artificially influence the number of clusters while having high accuracy and stability in 2D or 3D space. Theoretically, the increase in data dimensionality facilitates accurate characterization of the target reservoirs, but it also makes the computational effort grow exponentially and the accuracy of clustering may decrease due to data redundancy. Therefore, evaluating the optimal clustering dimensions, i.e., determining the number of clustering attributes, is an important factor affecting the clustering results. As the data dimension increases, the inter-group distance increases, while the tightness of the data within the group is decreasing. The best clustering result we want is that the inter-group distance is larger while the intra-group data tightness is also higher. Therefore, we use their ratios as evaluation indicator SC to achieve better evaluation of data dimensionality.
where Dis is the sum of the Euclidean distance between "exemplars". The joint distribution characteristics between the elastic parameters can be estimated with the help of the EM (expectation-maximization) algorithm, and the maximum probability density is considered as a representative point for characterizing a class of data, i.e., "exemplars". S_cat indicates the dispersion of the same class of data, and where c is the number of clusters, D is the entire dataset and C i represents the ith class. C i is the variance vector of the ith class of samples, (D) is the vector of variances for all samples. ‖ i.e., the Euclidean distance from the "exemplars", Similarly, ‖ (D)‖ is the L2-parametrization of the vector (D) . It can be directly observed that the indicators defined here help to evaluate the optimal clustering dimension and guide the selection of elastic attributes. Note that since we need to calculate the Euclidean distance between "exemplars", we need to normalize the data to eliminate the effect of the magnitude.

Numerical analysis
This section mainly contains two aspects to illustrate the effectiveness and accuracy of the proposed method: (1) validation of the rock-physics model and the sensitivity evaluation of elastic parameters; (2) verification of the feasibility of AP clustering strategy using logging data.
(16) SC= Dis S_cat Relative variation of elastic parameters when hydrate saturation changes from 10 to 50% Fig. 8 Analysis of the ability to distinguish hydrate reservoirs and water-bearing sediments by different dimensions of elastic attributes, red dots represent water-bearing sandstone, green dots represent hydrate layer. a 1D-shear modulus, b 2D-shear modulus and Young's modulus, c 3D-shear modulus, Young's modulus, and S-wave velocity, d 4D-shear modulus, Young's modulus, S-wave velocity, and S-impedance (the size of the point is denoted as the fourth spatial dimension of the S-impedance composition)

Rock-physics model validation and sensitivity evaluation of reservoir elastic parameters
Logging data from a hydrate work area in the South China Sea is used to verify the validity of the rock-physics model, and the specific logging curves are shown in Fig. 1. The core sampling proves that the hydrate reservoir is located at 195-220 mbsf. The apparent anomalies of high resistivity and velocity, as well as the small variations in lithology (GR), density, and porosity, are consistent with the typical characteristics of the presence of hydrates. Noted that the porosity curve is obtained from density logging and hydrate saturation is gained by the Archie equation calculation. The main elastic parameters applied in the modeling process are shown in Table 2.
The shear modulus of hydrates is considered in the modeling process, and we first analyze the response of rock elastic parameters to the shear modulus of hydrates. Figure 2 illustrates that the shear properties of hydrates do not affect the bulk modulus of saturated rocks, but the effect on the shear modulus is not negligible. This is because we regard the microscopic existence type of hydrate as the pore-filling mode, and therefore the properties of hydrates have directly influenced the characteristics of the equivalent fluid.
The consolidation parameter introduced in the modeling process has a great influence on the elastic parameters, especially in unconsolidated sediments. In addition, the estimation of consolidation parameter using the optimization method requires a reasonable initial value. Therefore, we combined the calculated values of Biot-Gassmann theory 5 Page 10 of 17 and the measured velocities of the logging data to set the initial values of the consolidation parameters by crossplot. Figure 3 indicates that the initial values of consolidation parameter for the hydrate target layer and the water-bearing layer are set to 20.0 and 25.0, respectively. Figure 4 shows the results of the P-and S-wave velocity estimation of hydrate reservoir using the proposed modeling method. The error analysis of velocity estimation (relative error) found that the error of P-and S-velocity prediction results are within 10% (Fig. 5a). Overall, the model predicts the velocities as required, and the estimation results of consolidation parameter obtained from the MC method indicate that the hydrate reservoir consolidation parameters are in the range of 19-24 (Fig. 5b), which is consistent with the previous understanding.
We analyzed the variation rule of elastic parameters with hydrate saturation (Fig. 6). Figure 7 shows the relative change in the elastic parameters as the hydrate saturation varies from 10 to 50%, from which the sensitivity ranking of the hydrate reservoir elastic parameters can be obtained. 15%-25% Sh <15% Fig. 10 Evaluating the effectiveness of AP clustering methods using P-and S-wave velocity data. a P-and S-wave velocity crossplot, b data input, c clustering results, d crossplot of hydrate saturation and velocity interpreted by logging (the arrow points to an increase in hydrate saturation) Next, the clustering dimension evaluation indicator SC is used to optimize the elastic attributes. As previously mentioned, the elastic parameters are first normalized to ensure that the effects of the magnitudes of the individual elastic parameters are eliminated in the calculation of the Euclidean distances. Secondly, the maximum probability distribution points of the data are used as "exemplars" to calculate the indicator SC . We tested the ability to distinguish hydrate reservoirs for 1D to 4D data consisting of shear modulus, Young's modulus, S-wave velocity and S-impedance attributes (Fig. 8), and the indicator SC is shown in Fig. 9. It can be found that the distance between "exemplars" increases as the data dimension rises. However, due to the increased dispersion of data, the growth rate of indicators is less than that of distance. When the data is 4-D, the evaluation indicator decreases obviously, which is mainly caused by the low sensitivity of S-impedance to the hydrate reservoir. Therefore, the data dimension is not simply the higher the better but depends on the sensitivity of the elastic attribute to the target reservoir. Ideally, high-dimensional data composed of elastic attributes with high sensitivity to the target reservoir would be preferred. In summary, we selected the 3D data space consisting of shear modulus, Young's modulus, and S-wave velocity to carry out accurate identification of hydrate reservoirs.

Evaluation of the effectiveness of AP clustering method on logging data
The data of well shA is used to test the effectiveness of AP clustering method for hydrate reservoir identification. The velocity of the hydrate reservoir is anomalous compared to the hydrocarbon-bearing sediments, which is the basis for the preliminary determination of hydrate presence using high-speed characteristics. Therefore, we first selected normalized P-and S-wave velocities as the data source to examine the clustering results. Figure 10a shows the crossplot of the P-and S-wave velocities. The 2D data composed of velocities are used as the original input data (Fig. 10b), and the result after clustering by AP strategy is shown in Fig. 10c. Since the main influencing factor of the velocity characteristics of hydrate reservoir is saturation, the clustering result is also in good agreement with the log interpretation of hydrate saturation when reflecting the abnormal velocity of hydrate reservoir (Fig. 10d). Similarly, we use predetermined shear modulus, Young's modulus, and S-wave velocity attributes combined with AP clustering to identify hydrate reservoirs. Figure 11 illustrates the clustering process and results, where the 3D data points are clearly divided into four clusters. To illustrate the stability and uncertainty of the clustering strategy, the number of repeated implementations shows that the clustering results remain almost consistent (Fig. 12). Meanwhile, the clustering results match well with the various characteristics of the logging curves (Fig. 13). For Fig. 13e, we matched the lithology with the clustering results by logging interpretation. Since the lithology of hydrate-bearing sediments is not changed (GR curve), the clustering results mainly reflect the differences in elastic response due to variations in reservoir physical parameters and pore fillings. The above test results show that the multiple elastic attribute fusion method based on AP clustering strategy can yield satisfactory hydrate reservoir identification results.

Field application
The proposed hydrate reservoir identification method is applied to the field area in the South China Sea. Since the Pliocene, seafloor slumping in the study area has been significant. The slippage fault inside the slump body provides a sufficient gas source for hydrate enrichment. The 3D seismic data and drilling results show that the sedimentary layers containing hydrates are located in the range of 150 to 250 mbsf. To meet the requirements of pre-stack seismic inversion, the CRP gathers are used as raw data, and partial angle stack seismic data are obtained by processing such as aberration removal, denoise, dip moveout, static correction, velocity analysis, and time-depth conversion. According to the depth of the reservoir and the quality of the seismic data, the incidence angles of partial angle stack seismic data are divided into 5°, 15°, and 25°, respectively. The predicted results of shear modulus, Young's modulus, and S-wave velocity are obtained by the inversion method (Appendix), while the accuracy of the elastic parameter predictions is confirmed by comparison with the logging curves (Fig. 14). The result of multiple elastic attribute clustering is shown in Fig. 15, and the hydrate reservoir identification result matches the hydrate development location confirmed by drilling. The comparison with the inversion results of the elastic parameters shows that the clustering strategy integrates the ability of each elastic attribute to best identify hydrate reservoirs, thus reducing the ambiguity of a single parameter for target reservoir identification. It should be pointed out that the inversion result data of elastic attributes do not have sufficient high-frequency information compared to the logging data, and therefore have less detail in the clustering results than in Fig. 13. In order to get a comprehensive understanding of the hydrate reservoir development in the whole work area, we mapped the slice along the horizon (Fig. 16). The multiple elastic attribute clustering results accurately identify the position of hydrate spatial development while providing better continuity (Fig. 16d).
In conclusion, the proposed multiple elastic attribute fusion method based on AP clustering strategy in hydrate reservoir Fig. 13 Comparison of AP clustering results and logging curves. a Resistivity, b shear modulus, c Young's modulus, d S-wave velocity, e clustering results (determine the corresponding relationship between clustering results and pore fluid characteristics according to logging interpretation) identification method is successfully applied, which illustrates its reasonableness and reliability.

Discussion
In this study, we implemented multiple elastic attribute clustering using AP clustering algorithm and obtained satisfactory results in hydrate reservoir identification. The high-dimensional data space composed of multiple elastic attributes is conducive to a comprehensive description of reservoir elastic properties, overcoming the deficiency of a single elastic attribute for hydrate reservoir characterization. The proposed multiple elastic attribute clustering method to identify hydrate reservoirs requires attention to several key aspects: (1) evaluation of reservoir-sensitive parameters. The reasonableness of the rock-physics model as a bridge to evaluate the sensitivity of elastic parameters plays a decisive role. In addition, direct rendezvous of logging data is a common method. Therefore, it is more reasonable to combine theoretical analysis and actual data to evaluate the sensitivity of elastic parameters. (2) Determine the optimal clustering dimension for the composition of elastic attributes. Although from a pure data perspective, the increase in data dimensionality facilitates a comprehensive description of the properties of the target reservoir, the sensitivity of different elasticity parameters has a variable weighting on the clustering results. Therefore, elastic attributes with more prominent sensitivity should be selected as the data source for clustering as much as possible. (3) The accurate inversion of the elastic parameters is extremely important. Due to the datadriven property of AP clustering strategy, the data source is the strongest influence on the clustering results. Therefore, inversion methods with better robustness and accuracy are aspects that deserve to be noticed.

Conclusions
We proposed a multiple elastic attribute fusion method based on AP clustering strategy to identify hydrate reservoirs. Firstly, the sensitivity of the elastic parameters was evaluated by a constructed rock-physics model considering the hydrate shear modulus. Secondly, the elastic attributes involved in the clustering were determined using the defined clustering best dimension indicator. Then, the accurate prediction results of the selected elastic properties were obtained in combination with the pre-stack seismic inversion method to provide the foundation for the data source. Finally, the fusion of different elastic attributes was achieved by clustering method, and multiple clusters reflecting the reservoir properties were obtained, which ultimately realized the accurate identification of hydrate reservoirs. The feasibility and validity of the proposed method are confirmed by field data. The method proposed in this paper also provides new ideas in identifying other types of reservoirs. In addition, the AP clustering strategy will be promising for geophysical fields such as seismic facies analysis, and engineering sweet spots forecast.

Bayesian pre-stack seismic inversion
Pre-stack inversion of elastic parameters is the data source for identification hydrate reservoir based on multiple elastic attribute clustering strategy. Bayesian pre-stack seismic inversion method containing smooth background constraints and boundary constraints is applied to elastic parameter prediction. The approximate equation for the P-wave reflectivity can be simplified as where N is the number of model parameters, m is model parameter, a i ( ) is weighting coefficient related to the angle of incidence. The posterior probability distribution p( m|S) of the model parameters in the Bayesian inversion framework is where N is the number of sampling points per seismic trace, m is the estimated model parameter, G is the wavelet matrix, 2 n and 2 k are the variance of the noise and model parameters, respectively. Take the maximum a posteriori probability of Eq. (19), yields Add smoothing background constraints to the objective function, where and i is the constraint factor of model parameters, × ln m i ∕m i r i = 1, 2, 3, the subscript r represents the model initial value. To avoid singular values in the inversion results, a boundary constraint is added to Eq. (21) using the Lagrange multiplier method.
where g i (m) is the boundary constraint function, is penalty parameter, is multiplicative vector. The parameter estimates can be obtained by iteratively solving the objective function. We give the iterative formulas and termination conditions as follows Accurate estimation of elastic parameters reflecting reservoir properties can be achieved using the above method.