AllAboutShale.com

Brittleness of rock plays a significant role in exploration and development of shale gas reservoirs. Young’s modulus and Poisson’s ratio are the key parameters for evaluating the rock brittleness in shale gas exploration because their combination relationship can quantitatively characterize the rock brittleness. The high-value anomaly of Young’s modulus and the low-value anomaly of Poisson’s ratio represent high brittleness of shale. The technique of pre-stack amplitude variation with angle inversion allows geoscientists to estimate Young’s modulus and Poisson’s ratio from seismic data. A model constrained basis pursuit inversion method is proposed for stably estimating Young’s modulus and Poisson’s ratio. Test results of synthetic gather data show that Young’s modulus and Poisson’s ratio can be estimated reasonably. With the novel method, the inverted Young’s modulus and Poisson’s ratio of real field data focus the layer boundaries better, which is helpful for us to evaluate the brittleness of shale gas reservoirs. The results of brittleness evaluation show a good agreement with the results of well interpretation.

Shale gas is a very important type of unconventional resource. The term refers to the unconventional gas stored in shale reservoirs. With the development of seismic exploration, a large amount of practice in unconventional shale reservoirs indicated that rock brittleness is one of the critical parameters to be taken into consideration in the evaluation of hydraulic fracturing. The study of shale brittleness is very important for shale gas exploration and development. An empirical brittleness cut-off is defined based on Young’s modulus and Poisson’s ratio (Grieser and Bray 2007; Rickman et al. 2008) as they control the relationship between stress and strain given by Hooke’s Law (Sena et al. 2011). A high-value anomaly of Young’s modulus and a low-value anomaly of Poisson’s ratio can be used to evaluate the rock brittleness and to infer ‘‘sweet spots’’ of shale gas reservoirs (Harris et al. 2011; Zong et al. 2013). Seismic inversion is the fundamental scientific tool used to obtain parameters concerning lithology and physical properties (Yin et al. 2015). Therefore, the estimation of Young’s modulus and Poisson’s ratio from pre-stack seismic data is a helpful guide for evaluating the brittleness of shale.

Amplitude variation with angle (AVA) inversion can be used to estimate the subsurface elastic properties from pre-stack seismic reflection data. However, the geophysical inversion problem in nature is an ill-conditioned problem because slight noise contained in the observed data will lead to enormous changes in the estimated parameters. Another problem of AVA inversion is that there are many models adequately fitting the data because the seismic data are band limited. It is common to add additional constraints to stabilize the inversion process and to reduce the number of solutions. This is generally referred to as regularization.

The regularization method was proposed by Tikhonov (1963), and the L-curve (Hansen 1992) was presented for selecting the regularization parameters which balance the data fitting term and trade-off function. A Bayesian approach is another method for stabilizing the inversion performance that treats the model parameter as a random variable with a probability distribution (Duijndam 1987; Buland and Omre 2003; Tarantola 2005; Yuan and Wang 2013; Yin and Zhang 2014). We seek the parameter esti- mates with a maximum posterior distribution combined with prior information of model parameters and the like-lihood function. The prior information can be the probability distribution of the model parameters and the geological information. In special cases, the regularization function is equivalent to the prior information in the Bayesian method.

The sparse solutions are full band; therefore, the sparse estimations are often viewed as high-resolution estimations (Levy and Fullagar 1981; Sacchi 1997; Alemie and Sacchi 2011). The inversion results via l2 norm regularization (Tikhonov 1963) or assumption of Gaussian probability distribution (Downton 2005) do not lead to high resolution because the estimates lack sparsity. The sparse reflection coefficients generate the blocky layer elastic parameters which suppress the side lobes. The development of the theory of sparse representation promotes the sparse inversion method. Theune et al. (2010) investigated the Cauchy and Laplace statistical distributions for their potential to recover sharp boundaries between adjacent layers. Based on the reflection dipole decomposition described by Chopra et al. (2006), Zhang et al. (2009, 2011) studied the basis pursuit inversion (BPI) of post and pre-stack seismic data, respectively, and got the sparse reflection coefficients and blocky layer elastic parameters, which is a high-resolution inversion method. Perez et al. (2013) proposed a hybrid Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) least-squares strategy that inverts the location of reflection by the FISTA algorithm (Beck and Teboulle 2009) first and then reevaluates the sparse (high resolution) reflection coefficients. However, this lacks the low-frequency information if we only use the seismic data and sparse regularization. The low-frequency information should be incorporated into the objective function to enhance the meaning of the inversion results and meanwhile promote the stability of the inversion implementation (Yin et al. 2008, 2014; Zong et al. 2012a; Yuan et al. 2015).

The ultimate goal of pre-stack inversion is to obtain the elastic information that can be used for evaluating hydrocarbon potential and the brittleness of the reservoir and to infer ideal drilling locations of ‘‘sweet spots’’ (Sena et al. 2011). Different linear approximations (Aki and Richards

1980; Gray et al. 1999; Russell et al. 2011) of the Zoeppritz equation introduced by Zoeppritz (1919) help us to directly estimate the elastic parameters (e.g., P-wave velocity, S-wave velocity, Lame´ parameters, bulk modulus, density) in which we are interested. The Young’s modulus and Poisson’s ratio can be calculated from the P-wave velocity, S-wave velocity, and density which can be inverted directly via Aki-Richards approximation. The density is difficult to invert, which will have a deleterious influence on the estimation of Young’s modulus. Parameters estimated indirectly will bring in more uncertainty in the inversion results (Zhang et al. 2009). In order to estimate the Young’s modulus (Y), Poisson’s ratio (P) and density (D) directly, Zong et al. (2012b) derived the linear approximation equation based on Young’s modulus, Pois- son’s ratio, and density and inverted the elastic parameters by Bayesian framework via Cauchy distribution as prior information. The approximation can be named as the YPD approximation. Zong et al. (2013) reformulated the elastic impedance equation in terms of Young’s modulus, Pois- son’s ratio and density based on the YPD approximation, and introduced a stable inversion method named elastic impedance varying with incident angle inversion with damping singular value decomposition (EVA-DSVD) inversion. In this study, we propose a model constrained basis pursuit inversion method to estimate the Young’s modulus, Poisson’s ratio, and density with the YPD approximation. The model constraint term is added into the objective function through a Bayesian framework. We also take a decorrelation of model parameters before inversion. The introduced model constraint promotes the stability of the inversion. Basis pursuit ensures the sparsity of reflection coefficients and the blocky structure of layer parameters. Model synthetic gather data with different signal-to- noise ratios are studied to test the proposed inversion method. The application on real data from shale reservoirs shows that Young’s modulus and Poisson’s ratio inverted by the proposed inversion method are reasonable for brittleness evaluation. The result of brittleness evaluation fits well with well interpretation.

Young’s modulus and Poisson’s ratio are the central parameters in predicting the brittleness of the subsurface layers. Young’s modulus can characterize the rigidity or brittleness of rocks, and Poisson’s ratio can be regarded as a kind of fluid factor which can be used for pore fluid identification. The reflection coefficients approximate equation was derived in terms of Young’s modulus, Poisson’s ratio, and density (YPD approximation) with the hypothesis of planar incident wave by Zong et al. (2012c):

where θ is the incident angle; R(θ) is the reflection coefficients; k stands for the square of the average S-to-P velocity ratio; and ΔE/E, Δr/r, and Δq/q represent the reflection coefficients of Young’s modulus, Poisson’s ratio, and density, respectively.

In order to perform the inversion for Young’s modulus, Poisson’s ratio and density, we should firstly build the forward model. Combining the convolution model and YPD approximation, we can get the pre-stack data in angle domain shown as Eq. (2).

where W(θi) represents the wavelet matrix and CE (θi), Cθ (θi), and Cq (θi) represent the weighting coefficients of Young’s modulus reflectivity vector rE, Poisson’s ratio reflectivity vector rθ, and density reflectivity vector rρ, respectively. The product of the wavelet matrix and weighting coefficients makes up the kernel matrix G. Set-equation can be written in a linear equation as d = Gr.

Decorrelation of model space parameters can enhance the stability of the three-parameter AVA inversion (Downton2005; Zong et al. 2012b). We took the singular value decomposition (SVD) for covariance matrix Cr of model elastic parameters:

where θ2 is the variance of Young’s modulus; θEθ is covariance of Young’s modulus and Poisson’s ratio, and so on; v is a matrix made up of three eigenvectors; and S is the diagonal matrix made up of the positive decreasing singular values.

For N samples, the inverse of the decorrelation matrix V can be expressed as the Kronecker product of v and an N-order identity matrix. Therefore, the inverse of decorrelation matrix V is expressed as Eq. (4):

In this case, the kernel matrix G becomes G=GV-1, and the model vector becomes r= Vr. Therefore, the forward modeling equation should be written as

Then we used the dipole reflection coefficient decomposition method to update the forward model of pre-stack BPI inversion for Young’s modulus, Poisson’s ratio, and density. The reflection coefficient decomposition method is shown in Fig. 1. In this case, the vector of reflection coefficients containing Young’s modulus, Poisson’s ratio, and density can be written as Eq. (6).

where D stands for the dipole reﬂectivity decomposition operator. mE, mθ and mρ are the sparse coefﬁcients of Young’s modulus, Poisson’s ratio, and density, respectively, under the reﬂectivity decomposition, mE=[aET bET], mθ=[aθT bθT] and mρ=[aρT bρT]. D is a large scale matrix consisting of three reﬂectivity decomposition operators. m is the sparse coefﬁcients vector consisting of mE, mr, and mρ. We can obtain the forward model shown as Eq. (7) by substituting Eq. (6) into Eq. (5):

In this study, we constructed the objective function of the AVA inversion under the Bayesian framework (Ulrych et al. 2001). The Bayesian theorem is given by

where p (m/d) is the posterior probability distribution function (PDF), p (d/m) is the likelihood function that is the PDF of noise, p (m) is the prior information of the parameter m, and p (d) is the normalization factor which can be ignored as it is a constant value. Therefore, the Bayesian theorem is often expressed as Eq. (9) without the scaling factor:

We suppose that the noise obeys a Gaussian distribution; hence, the likelihood function should be

where Xd is the noise covariance matrix. For simplicity, we suppose that the noise is uncorrelated, so the covariance matrix should be Xd=θd2Id , where Id is the identity matrix and θd2 is the variance of the Gaussian distributed noise.

In this paper, we assume that the prior distribution is constructed by two terms:

where the first term, pt (m), is the probability distribution of m which represents the sparsity of the coefficients and the second term, pl (m), is the low-frequency model information that can enhance the lateral continuity. The prior information stabilized the inversion process and provided a principle to choose the ‘‘best’’ solution that can adequately fit the observed data.

In order to recover the discontinuity of the layer properties, the minimum l1 norm that works well in selecting the sparse solution practically should be taken into consideration to constrain the inversion. This l1 norm regularization can be incorporated into the Bayesian approach as the Laplacian distribution with a mean of zero:

In the lateral term, we suppose that the error between the inversion and low-frequency model obeys a normal distribution.

where X(m,ξ) is the covariance matrix associated with the three elastic properties: Young’s modulus, Poisson’s ratio, and density. Here we assume that the parameters at each sample are independent as we took decorrelation of the models parameters, and then X(m,ξ)= θ2(m,ξ)I(m,ξ) where θ2(m,ξ) is the variance of the coefﬁcients for estimation, and I(m,ξ) is the identity matrix. ξ is the vector made up of the relevant Young’s modulus, Poisson’s ratio, and density; C is made up of diagonal integrated matrix

The expression of ξ and C can be written as Eq. (14):

Under Bayes’ framework, we can estimate the solution as the maximum a posteriori (MAP) solution. Substituting Eq. (12) and Eq. (13) into Eq. (11), the prior information can be written as

Substituting the likelihood function Eq.(10) and prior distribution Eq. (15) into the Bayesian theorem Eq. (9), we get the objective function shown as Eq. (16) under the Bayesian inference:

θ where, λ and μ are the trade-off factors which balance the overall impact of the regularization, and λ=θd2/θm λ=θd2/θ2(m,ξ).

The objective function Eq. (16) can be viewed as the normal expression for a basis pursuit problem via an augmented matrix:

Accordingly, utilizing the Gradient Projection for Sparse Reconstruction (GPSR) method (Figueiredo et al. 2007), we minimized the objective function J(m) to obtain the sparse estimates. After that, the reflection coefficients with isolated spikes can be obtained by Eq. (6). The inverted results contain low frequencies as the low-frequency trend model data were added into the objective function as a distribution.

We tested the validity of our proposed inversion method with well log data. The angle gather data with free noise (Fig. 2a) were synthesized with Zoeppritz equations in the time domain by utilizing the real well logs of P-wave velocity, S-wave velocity, and density and a 35 Hz Ricker wavelet for incident angles ranging from 0° to 35°. Figure 2b–d displays the log curves of Young’s modulus, Poisson’s ratio, and density. The blue curves shown in Fig. 2b–d are the real models and the red curves are the inversion results. From Fig. 2b–d, we can clearly see that the Young’s modulus, Poisson’s ratio, and density can be inverted reasonably with free noise. The error of the inverted density is a little bigger than that of the other two parameters.

Fig. 2 a The synthetic seismic angle gather with free noise; b The estimated Young’s modulus; c The estimated Poisson’s ratio; d The estimated density. b–d the red lines mean the inverted results and the blue lines represent the real model curves

In order to verify the stability of the inversion method, we added random Gaussian noise to the synthetic gather data with different signal-to-noise ratios (SNRs). The gather traces are displayed in Figs. 3a, 4a and 5a, and the SNRs are 4:1, 2:1, and 1:1, respectively. The inversion results in different gather traces are shown in Figs. 3b–d, 4b–d and 5b–d, respectively.

Fig. 3 a The synthetic seismic angle gather with SNR=4:1; b The estimated Young’s modulus; c The estimated Poisson’s ratio; d The estimated density. b–d the red lines mean the inverted results, and the blue lines represent the real model curves

It is very clear that the inversion results estimated by the proposed inversion method match well with the real models as the low-frequency model enables the inversion results to approximate the real models. Although the inversion results are influenced by the noise, especially the density, the Young’s modulus and Poisson’s ratio can match well with the real models to some degree so that the inversion results are helpful for us to evaluate the brittleness of the layer.

Fig. 4 a The synthetic seismic angle gather with SNR = 2:1; b The estimated Young’s modulus; c The estimated Poisson’s ratio; d The estimated density. b–d the red lines mean the inverted results, and the blue lines represent the real model curves

Fig. 5 a The synthetic seismic angle gather with SNR = 1:1; b The estimated Young’s modulus; c The estimated Poisson’s ratio; d The estimated density. b–d the red lines mean the inverted results, and the blue lines represent the real model curves

The inversion method was applied to real partial angle stack seismic data, and the sampling interval of the seismic data is 2 ms. Figure 6 displays the used three partial angle stack seismic sections, and the mean angles of the seismic data in Fig. 6a–c are 8 , 16 , and 24 , respectively. In Fig. 6a–c, the green ellipse circles the target reservoir, and a well is drilled through the target at CDP 156. Figure 6d–f displays the initial models of Young’s modulus, Poisson’s ratio, and density, respectively. The initial models are established by spatial interpolation and extrapolation and low-pass filtering. The inverted isolated reflectivity spikes of rE, rσ, and rρ are shown in Fig. 7a–c, respectively. The structure of the reflectivity estimates is similar to the partial angle stack seismic profile and appears to have a better resolving power for the layer boundaries.

Fig. 6 The partial stack seismic proﬁle and initial model of Young’s modulus, Poisson’s ratio, and density. a Partial stack seismic data with a small incident angle range, the mean angle is 8°; b Partial stack seismic data with a medium incident angle range, the mean angle is 16°; c Partial stack seismic data with a large incident angle range, the mean angle is 24°; d–f The proﬁles of initial model of Young’s modulus, Poisson’s ratio, and density, respectively

The blocky Young’s modulus, Poisson’s ratio, and density displayed in Fig. 7d–f are obtained from the estimated sparse reflectivity by using Eq. (18). The blocky results from model constrained BPI focus on the layer boundaries well, which is useful for us to interpret the inversion results. To some extent, the lateral continuity is improved because the low- frequency trend model is continuous laterally. The calculated logs of Young’s modulus, Poisson’s ratio and density are inserted into the sections.

Fig. 7 The inverted reﬂection properties of Young’s modulus, Poisson’s ratio, and density and the properties generated by the reﬂections using Eq. (18). a Reﬂection of Young’s modulus, b Reﬂection of Poisson’s ratio, c Reﬂection of density, d Young’s modulus, e Poisson’s ratio, f Density

Figure 8a–c plot the inverted Young’s modulus, Poisson’s ratio, and density (red lines) at the well location, aligning with the well logs (dark lines). From Fig. 8, we can draw a conclusion that the inverted results of Young’s modulus and Poisson’s ratio have a good fit with the logs, while the density inversion does not match as well as the other two elastic parameters as the maximum angle is not large enough for us to invert the density information. From the inverted results shown in Fig. 7d–f, we can clearly see that the Young’s modulus exhibits high anomalous values and Poisson’s ratio shows low anomalous values in the target circled with dark ellipses, which means that the brittleness of the circled target is high. The drilling shows the high brittleness at 2.3 s. Therefore, the brittleness evaluated by the inverted Young’s modulus and Poisson’s ratio using the proposed inversion method is consistent with the drilling.

Fig. 8 The comparison between inversion results and well logs. a Young’s modulus, b Poisson’s ratio, c density

In this paper, we presented a novel stable inversion method to estimate the Young’s modulus and Poisson’s ratio for brittleness evaluation from pre-stack seismic data with the YPD approximation. We introduced the low-frequency trend model into the basis pursuit inversion implementation. Therefore, we can call this method model constrained BPI. We derived the objective function of model constrained BPI by the Bayesian theorem. In the improved method, the low- frequency trend model as prior knowledge stabilized the inversion and improved the lateral continuity because the low-frequency trend model is continuous in the space axis. The l1 norm of the model parameters and the GPSR algorithm kept the sparsity of inversion results so that we can obtain the isolated reflectivities and the elastic parameters with discontinuous jumps in the time axis. The model test results showed that we can obtain the high-precision Young’s modulus, Poisson’s ratio, and density. The real data application is performed to confirm the validity of the proposed method, and it showed that the high anomalous value of Young’s modulus and low anomalous value of Poisson’s ratio which mean high brittleness matched well with the brittleness interpretation of drilling.

We would like to acknowledge the sponsorship of the National ‘‘973 Program’’ of China (2013CB228604) and the National Grand Project for Science and Technology (2011ZX05030-004-002), China Postdoctoral Science Foundation (2014M550379), Natural Science Foundation of Shandong (2014BSE28009), Science Foundation for Post-doctoral Scientists of Shandong (201401018), and Science Foundation for Post-doctoral Scientists of Qingdao and Science Foundation from SINOPEC Key Laboratory of Geophysics (33550006-14-FW2099-0038). We also acknowledge the support of the Australian and Western Australian governments and the North West Shelf Joint Venture partners, as well as the Western Australian Energy Research Alliance (WA:ERA).

Aki K, Richards PG. Quantitative seismology. San Francisco: W. H. Freeman; 1980.

Alemie W, Sacchi MD. High-resolution three-term AVO inversion by means of a trivariate cauchy probability distribution. Geo- physics. 2011;76(3):R43–55.

Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imaging Sci.2009;2(1):183–202.

Buland A, Omre H. Bayesian linearized AVO inversion. Geophysics. 2003;68(1):185–98.

Chopra S, Castagna J, Portniaguine O. Seismic resolution and thin- bed reflectivity inversion. CSEG Rec. 2006;31(1):19–25.

Downton JE. Seismic parameter estimation from AVO inversion. Ph.D. Thesis, University of Calgary. 2005.

Duijndam AJW. Detailed Bayesian inversion of seismic data. Ph.D. Thesis, Technische Universiteit Delft. 1987.

Figueiredo MAT, Nowak RD, Wright SJ. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems. IEEE J Sel Top Signal Process. 2007;1(4):586–97.

Gray D, Goodway B, Chen T. Bridging the gap: using AVO to detect changes in fundamental elastic constants. SEG Technical Program Expanded Abstracts. 1999.

Grieser WV, Bray JM. Identification of production potential in unconventional reservoirs. Production and Operations Symposium. 2007.

Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. Soc Ind Appl Math. 1992;34(4):561–80.

Harris NB, Miskimins JL, Mnich CA. Mechanical anisotropy in the Woodford Shale, Permian Basin: origin, magnitude, and scale. Lead Edge. 2011;30(3):284–91.

Levy S, Fullagar PK. Reconstruction of a sparse spike train from a portion of its spectrum and application to high-resolution deconvolution. Geophysics. 1981;46(9):1235–43.

Perez DO, Velis DR, Sacchi MD. High-resolution prestack seismic inversion using a hybrid FISTA least-squares strategy. Geophysics. 2013;78(5):R185–95.

Rickman R, Mullen MJ, Petre J E, et al. A practical use of shale petrophysics for stimulation design optimization: All shale plays are not clones of the Barnett shale. SPE Annual Technical Conference and Exhibition. 2008.

Russell BH, Gray D, Hampson DP. Linearized AVO and poroelas- ticity. Geophysics. 2011;76(3):C19–29.

Sacchi MD. Reweighting strategies in seismic deconvolution. Geo- phys J Int. 1997;129(3):651–6.

Sena A, Castillo G, Chesser K, et al. Seismic reservoir characteri zation in resource shale plays: ‘‘Sweet spot’’ discrimination and optimization of horizontal well placement. SEG Technical Program Expanded Abstracts. 2011.

Tarantola A. Inverse problem theory and methods for model parameter estimation. Soc Ind Appl Math. 2005;9(5):1597–1620. doi:10.1137/1.9780898717921.

Theune U, Jensa˚s IØ, Eidsvik J. Analysis of prior models for a blocky inversion of seismic AVA data. Geophysics. 2010;75(3): C25–35.

Tikhonov AN. Solution of incorrectly formulated problems and the regularization method. Sov Math Dokl. 1963;4:1035–8.

Ulrych TJ, Sacchi MD, Woodbury A. A Bayes tour of inversion: a tutorial. Geophysics. 2001;66(1):55–69.

Yin XY, Zhang SX. Bayesian inversion for effective pore-fluid bulk modulus based on fluid-matrix decoupled amplitude variation with offset approximation. Geophysics. 2014;79(5):R221–32.

Yin XY, Yang PJ, Zhang GZ. A novel prestack AVO inversion and its application. SEG Technical Program Expanded Abstracts. 2008.

Yin XY, Zong ZY, Wu GC. Seismic wave scattering inversion for fluid factor of heterogeneous media. Sci China Earth Sci. 2014;57(3):542–9.

Yin XY, Zong ZY, Wu GC. Research on seismic fluid identification driven by rock physics. Sci China Earth Sci. 2015;58(2):159–71.

Yuan SY, Wang SX. Spectral sparse Bayesian learning reflectivity inversion. Geophys Prospect. 2013;61(4):735–46.

Yuan SY, Wang SX, Luo CM, et al. Simultaneous multitrace impedance inversion with transform-domain sparsity promotion. Geophysics. 2015;80(2):R71–80.

Zhang R, Castagna J. Seismic sparse-layer reflectivity inversion using basis pursuit decomposition. Geophysics. 2011;76(6):R147–58. Zhang SX, Yin XY, Zhang FC. Fluid discrimination study from fluid elastic impedance (FEI). SEG Technical Program Expanded bstracts. 2009.

Zhang R, Sen MK, Srinivasan S. A prestack basis pursuit seismic inversion. Geophysics. 2013;78(1):R1–11.

Zoeppritz K. VII b. U¨ ber Reflexion und Durchgang seismischer

Wellen durch Unstetigkeitsfla¨chen. Nachrichten von der Koniglichen Gesellschaft der Wissenschaften zu Go¨ ttingen, Mathematisch-physikalische Klasse. 1919; 66–84.

Zong ZY, Yin XY, Wu GC. AVO inversion and poroelasticity with P- and S-wave moduli. Geophysics. 2012a;77(6):N17–24.

Zong ZY, Yin XY, Wu GC. Elastic impedance variation with angle inversion for elastic parameters. J Geophys Eng. 2012b;9(3): 247–60.

Zong ZY, Yin XY, Zhang F, et al. Reflection coefficient equation and pre-stack seismic inversion with Young’s modulus and Poisson ratio. Chin J Geophys. 2012c;55(11):3786–94 (in Chinese).

Zong ZY, Yin XY, Wu GC. Elastic impedance parameterization and inversion with Young’s modulus and Poisson’s ratio. Geo- physics. 2013;78(6):N35–42.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creative commons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Author: Xing-Yao Yin xyyin@upc.edu.cn