- Home
- A-Z Publications
- Geophysical Prospecting
- Previous Issues
- Volume 67, Issue 7, 2019
Geophysical Prospecting - Volume 67, Issue 7, 2019
Volume 67, Issue 7, 2019
-
-
Time‐frequency decomposition of seismic signals via quantum swarm evolutionary matching pursuit
More LessABSTRACTMatching pursuit belongs to the category of spectral decomposition approaches that use a pre‐defined discrete wavelet dictionary in order to decompose a signal adaptively. Although disengaged from windowing issues, matching point demands high computational costs as extraction of all local structure of signal requires a large size dictionary. Thus in order to find the best match wavelet, it is required to search the whole space. To reduce the computational cost of greedy matching pursuit, two artificial intelligence methods, (1) quantum inspired evolutionary algorithm and (2) particle swarm optimization, are introduced for two successive steps: (a) initial estimation and (b) optimization of wavelet parameters. We call this algorithm quantum swarm evolutionary matching pursuit. Quantum swarm evolutionary matching pursuit starts with a small colony of population at which each individual, is potentially a transformed form of a time‐frequency atom. To attain maximum pursuit of the potential candidate wavelets with the residual, the colony members are adjusted in an evolutionary way. In addition, the quantum computing concepts such as quantum bit, quantum gate, and superposition of states are introduced into the method. The algorithm parameters such as social and cognitive learning factors, population size and global migration period are optimized using seismic signals. In applying matching pursuit to geophysical data, typically complex trace attributes are used for initial estimation of wavelet parameters, however, in this study it was shown that using complex trace attributes are sensitive to noisy data and would have lower rate of convergence.
The algorithm performance over noisy signals, using non‐orthogonal dictionaries are investigated and compared with other methods such as orthogonal matching pursuit. The results illustrate that quantum swarm evolutionary matching pursuit has the least sensitivity to noise and higher rate of convergence. Finally, the algorithm is applied to both modelled seismograms and real data for detection of low frequency anomalies to validate the findings.
-
-
-
3D Seismic survey design using mixed‐radix based algorithm inversion
Authors Atahebson B. Santos, Milton J. Porsani and Manuelle S. GóisABSTRACTThe determination of three‐dimensional geometry and acquisition parameters, the seismic acquisition survey design, is constantly subject of studies in obtaining data with the highest seismic quality, operational efficiency and cost minimization. In this paper, we propose a methodology for inverting geometry parameters of three‐dimensional orthogonal land seismic surveys based on a direct search method using a mixed‐radix based algorithm. In this algorithm, the search space is discretized on a mixed‐radix base, which depends on the extreme values and the search resolution of each parameter. We will show how to reparametrize the orthogonal acquisition geometry elements in order to obtain the independents and integers parameters that are necessary to construct the mixed‐radix base. For the optimization purpose, we define an objective function to contemplate target parameters associated with the elements of the acquisition geometry directly related to the geophysical and operational constraints. Taking in account that the mathematical functions and the objective function we define for the problem have no significant computational cost, all model space parameters are fast and efficiently tested. We applied the algorithm, using as input data, provided by a one‐line roll orthogonal reference geometry, assuming a pair of geological objectives as shallow and deep targets. All selected models that meet both the proposed objectives and the constraints are organized by decreasing order of fitness so that with the mixed‐radix inversion algorithm we found not only the best model, but also a set of suitable models. Likewise, with the best set of geometries, it is possible to establish a direct comparison between them, analysing their adherence to the technical and operational requirements according to the availability and degree of detail of each one. We show the top 10 best results as a table, allowing a direct comparison between all aspects of these geometries, and we summarize the results showing graphically the fitness of all selected geometries and the inverted geometry elements for the 1000 best geometries. These graphical displays provide a direct way to understand how each model behaves as the fitness decreases. The algorithm is very flexible and its application can be extended to any environment and type of acquisition geometry, and in any phase study of an area be it regional, exploratory or development.
-
-
-
Perfectly matched layer boundary conditions for frequency‐domain acoustic wave simulation in the mesh‐free discretization
Authors Xin Liu, Yang Liu, Zhiming Ren and Bei LiABSTRACTMesh‐free discretization, flexibly distributing nodes without computationally expensive meshing process, is able to deal with staircase problem, oversampling and undersampling problems and saves plenty of nodes through distributing nodes suitably with respect to irregular boundaries and model parameters. However, the time‐domain mesh‐free discretization usually exhibits poorer stability than that in regular grid discretization. In order to reach unconditional stability and easy implementation in parallel computing, we develop the frequency‐domain finite‐difference method in a mesh‐free discretization, incorporated with two perfectly matched layer boundary conditions. Furthermore, to maintain the flexibility of mesh‐free discretization, the nodes are still irregularly distributed in the absorbing zone, which complicates the situation of artificial boundary reflections. In this paper, we implement frequency‐domain acoustic wave modelling in a mesh‐free system. First, we present the perfectly matched layer boundary condition to suppress spurious reflections. Moreover, we develop the complex frequency shifted–perfectly matched layer boundary condition to improve the attenuation of grazing waves. In addition, we employ the radial‐basis‐function‐generated finite difference method in the mesh‐free discretization to calculate spatial derivatives. The numerical experiment on a rectangle homogeneous model shows the effectiveness of the perfectly matched layer boundary condition and the complex frequency shifted–perfectly matched layer boundary condition, and the latter one is better than the former one when absorbing large angle incident waves. The experiment on the Marmousi model suggests that the complex frequency shifted–perfectly matched layer boundary condition works well for complicated models.
-
-
-
Sparse Bayesian linearized amplitude‐versus‐angle inversion
Authors Amir Mollajan, Hossein Memarian and Beatriz QuintalABSTRACTThe technique of seismic amplitude‐versus‐angle inversion has been widely used to estimate lithology and fluid properties in seismic exploration. The amplitude‐versus‐angle inversion problem is intrinsically ill‐posed and generally stabilized by the use of L2‐norm regularization methods but with drawback of smoothing important boundaries between adjacent layers. In this study, we propose a sparse Bayesian linearized solution for amplitude‐versus‐angle inversion problem to preserve the sharp geological interfaces. In this regard, a priori constraint term with two regularization functions is presented: the sparse constraint regularization and the low‐frequency model information. In addition, to obtain high‐resolution reflectivity estimation, the model parameters decorrelation technique combined with dipole decomposition method is employed. We validate the applicability of the presented method by both synthetic and real seismic data from the Gulf of Mexico. The accuracy improvement of the presented method is also confirmed by comparing the results with the commonly used Bayesian linearized amplitude‐versus‐angle inversion.
-
-
-
Combined inversion of first‐arrival travel times and reflection travel times
Authors Chongjin Zhao, Luolei Zhang, Peng Yu, Yuzhu Liu and Shaokong FengABSTRACTIn this work, we propose a method for determining reflection travel times based on the acquisition of first‐arrival travel times via the fast sweeping method. The accuracy of this scheme was proven by conducting model experiments to establish a foundation for first‐arrival tomography, reflection tomography and combined tomography. Reflection tomography was subsequently achieved using the adjoint‐state method; on this basis, we propose a combined tomography method involving both first‐arrival and reflection tomography. In the model experiments, excellent results were obtained via first‐arrival tomography, reflection tomography and our combined tomography method. Finally, full‐waveform inversion was performed, with the inversion produced by combined tomography used as the initial model. Excellent results were obtained using this approach. However, combined tomography reproduced and characterized the model much more realistically.
-
-
-
Multiwell study of seismic attenuation at the CO2CRC Otway project geosequestration site: Comparison of amplitude decay, centroid frequency shift and 1D waveform inversion methods
Authors A. Pirogova, R. Pevzner, B. Gurevich, S. Glubokovskikh and K. TertyshnikovABSTRACTAt the CO2CRC Otway geosequestration site, the abundance of borehole seismic and logging data provides a unique opportunity to compare techniques of Q (measure of attenuation) estimation and validate their reliability. Specifically, we test conventional time‐domain amplitude decay and spectral‐domain centroid frequency shift methods versus the 1D waveform inversion constrained by well logs on a set of zero‐offset vertical seismic profiles. The amplitude decay and centroid frequency shift methods of Q estimation assume that a seismic pulse propagates in a homogeneous medium and ignore the interference of the propagating wave with short‐period multiples. The waveform inversion explicitly models multiple scattering and interference on a stack of thin layers using high‐resolution data from sonic and density logs. This allows for stable Q estimation in small depth windows (in this study, 150 m), and separation of the frequency‐dependent layer‐induced scattering from intrinsic absorption. Besides, the inversion takes into account band‐limited nature of seismic data, and thus, it is less dependent on the operating frequency bandwidth than on the other methods. However, all considered methods of Q estimation are unreliable in the intervals where subsurface significantly deviates from 1D geometry. At the Otway site, the attenuation estimates are distorted by sub‐vertical faults close to the boreholes. Analysis of repeated vertical seismic profiles reveals that 15 kt injection of the CO2‐rich fluid into a thin saline aquifer at 1.5 km depth does not induce detectable absorption of P‐waves at generated frequencies 5–150 Hz, most likely because the CO2 plume in the monitoring well is thin, <15 m. At the Otway research site, strong attenuation Q ≈ 30–50 is observed only in shaly formations (Skull Creek Mudstone, Belfast Mudstone). Layer‐induced scattering attenuation is negligible except for a few intervals, namely 500–650 m from the surface, and near the injection interval, at around 1400–1550 m, where Qscat ≈ 50–65.
-
-
-
Pattern‐guided dip estimation with plane‐wave destruction filters
Authors Zhiguang Xue, Houzhu Zhang, Yang Zhao and Sergey FomelABSTRACTWe propose to use pattern‐guided dip estimation to mitigate aliasing problem that possibly exists in structure‐oriented data processing. A straightforward and effective approach of generating pattern‐guided dip is presented, which generally involves three rounds of standard dip estimation with plane‐wave destruction filters. The first use of plane‐wave destruction filter is for generating a mask operator distinguishing aliased and non‐aliased data, based on measuring the uncertainty of the first dip estimation. The second plane‐wave destruction filter uses the aliasing‐free portions of the input data, and the dip in the aliasing‐affected area is automatically padded with the ‘pattern’ dip by smoothing regularization. The result of the second plane‐wave destruction filter is used as the initial dip for the inversion of the last‐pass plane‐wave destruction filter, which produces a pattern‐guided dip. For some specific applications, the mask operator can be easily generated through other methods, and we can skip the first dip estimation. Two numerical examples, related to picking information using predictive painting and structure‐oriented interpolation, respectively, demonstrate that our pattern‐guided dip can effectively mitigate the aliasing problem in structure‐oriented data processing.
-
-
-
A fast and robust two‐point ray tracing method in layered media with constant or linearly varying layer velocity
Authors Xinding Fang and Xiaofei ChenABSTRACTA fast and robust method for two‐point ray tracing in one‐dimensional layered media is presented. This method is applicable to layered models with constant or linearly varying isotropic layer velocity. For given model properties and source and receiver positions, a ray path can be uniquely determined once its ray parameter (i.e. horizontal slowness) is known. The ray parameter can be obtained by numerically solving the nonlinear offset (i.e. source–receiver horizontal distance) equation using Newton's method, which generally works well at near and mid offsets. However, Newton's method becomes hard to converge at large offsets due to the oversensitivity of offset to ray parameter. Based on the analysis of the characteristic of the offset equation, a modified ray parameter is proposed and used to replace the generic ray parameter in numerical calculation. Numerical experiments show that the iteration process becomes stable and converges rapidly with the modified ray parameter. Moreover, a rational function that asymptotically approximates the shape of the offset equation is introduced for obtaining good initial estimates of the modified ray parameter. Numerical tests show that this method is robust in any situation, and an accurate ray parameter can be obtained within two or three iterations for a wide range of model velocity structure and source–receiver distance. Furthermore, the proposed two‐point ray tracing method is easy to implement.
-
-
-
Mitigating the effects of sand dunes on seismic data from the Rub al Khali basin, Saudi Arabia
ABSTRACTThe Rub Al Khali region in Saudi Arabia is characterized by the presence of sand dunes separated by salt flats, also called Sabkhas. In general, the elevation of dunes in this region varies between 90 and 250 m above sea level. The presence of these sand dunes, along with the rapidly changing surface topography poses challenges for seismic data acquisition and processing. The high contrast in acoustic impedance between the dune base and the underlying formation often results in amplification of seismic waves that are recorded at stations located on the surface of sand dunes. Attempts to address the issue using conventional surface‐consistent amplitude scaling methods without reducing these amplification effects generally fail, thus compromising the suitability of the processed data for amplitude versus offset analysis. In this study, we propose a new reference site technique to reduce the effects of sand dune amplification, enabling the production of data sets that are suitable for amplitude versus offset processing. The proposed technique uses a deterministic approach to derive surface‐consistent, frequency‐dependent de‐amplification functions for shots and stations located on the dunes only. Two‐dimensional synthetic and field data examples show that the technique significantly reduces the effects of sand dune amplification.
-
-
-
Time‐frequency analysis based on curvelet transforms with time skewing
Authors Jiao Wang, Zhenchun Li, Lutz Gross and Stephen TysonABSTRACTOil and gas exploration gradually changes to the deep and complex areas. The quality of seismic data restricts the effective application of conventional time‐frequency analysis technology, especially in the case of low signal‐to‐noise ratio. To address this problem, we propose a curvelet‐based time‐frequency analysis method, which is suitable for seismic data, and takes into account the lateral variation of seismic data. We first construct a kind of curvelet adapted to seismic data. By adjusting the rotation mode of the curvelet in the form of time skewing, the scale parameter can be directly related to the frequency of the seismic data. Therefore, the curvelet coefficients at different scales can reflect the time‐frequency information of the seismic data. Then, the curvelet coefficients, which represent the dominant azimuthal pattern, are converted to the time‐frequency domain. Since the curvelet transform is a kind of sparse representation for the signal, the screening process of the dominant coefficient masks most of the random noise, which enables the method to adapt for the low signal‐to‐noise ratio data. Results of synthetic and field data experiments using the proposed method demonstrate that it is a good approach to identify weak signals from strong noise in the time‐frequency domain.
-
-
-
Research note: deblended‐data reconstruction using generalized blending and deblending models
ABSTRACTWe introduce a concept of generalized blending and deblending, develop its models and accordingly establish a method of deblended‐data reconstruction using these models. The generalized models can handle real situations by including random encoding into the generalized operators both in the space and time domain, and both at the source and receiver side. We consider an iterative optimization scheme using a closed‐loop approach with the generalized blending and deblending models, in which the former works for the forward modelling and the latter for the inverse modelling in the closed loop. We applied our method to existing real data acquired in Abu Dhabi. The results show that our method succeeded to fully reconstruct deblended data even from the fully generalized, thus quite complicated blended data. We discuss the complexity of blending properties on the deblending performance. In addition, we discuss the applicability to time‐lapse seismic monitoring as it ensures high repeatability of the surveys. Conclusively, we should acquire blended data and reconstruct deblended data without serious problems but with the benefit of blended acquisition.
-
-
-
Rock elasticity as a function of the uniaxial stress: laboratory measurements and theoretical modelling of vertical transversely isotropic and orthorhombic shales
Authors V.A. Sviridov, S.I. Mayr and S.A. ShapiroABSTRACTWe apply a rock‐physics model that describes the relationship between the effective stress and rock elasticity. We experimentally obtain and analyse a data set containing one vertical transversely isotropic and one orthorhombic shale sample. The vertical transversely isotropic symmetry of the first sample is caused by the layered structure of the rock. The seismic orthorhombicity of the second sample could be explained after microscopic analysis of thin section, which demonstrates an imperfect disorder of inhomogeneities. Both samples were loaded uniaxially in a quasi‐static regime. During the loading, we measured stress‐dependent seismic velocities and sample deformations. For the analysis of the stress‐dependent velocities and stiffnesses, we modelled the measured data set using a recent generalization of the porosity deformation approach. Comparison of the experimentally determined and numerically modelled data supports the applicability of the theory and helps in the interpretation of experimentally obtained data. In agreement with the theory, uniaxial stress increases the elliptic component of the seismic anisotropy and does not impact the anellipticity parameter. We demonstrate the distinct influence of the stiff and compliant porosities on the stress sensitivity of the elastic properties.
-
-
-
Linearized amplitude variation with offset and azimuth and anisotropic poroelasticity
Authors Xinpeng Pan, Guangzhi Zhang and Xingyao YinABSTRACTThe normal‐to‐shear weakness ratio is commonly used as a fracture fluid indicator, but it depends not only on the fluid types but also on the fracture intensity and internal architecture. Amplitude variation with offset and azimuth is commonly used to perform the fluid identification and fracture characterization in fractured porous rocks. We demonstrate a direct inversion approach to utilize the observable azimuthal data to estimate the decoupled fluid (fluid/porosity term) and fracture (normal and shear weaknesses) parameters instead of the calculation of normal‐to‐shear weakness ratio to help reduce the uncertainties in fracture characterization and fluid identification of a gas‐saturated porous medium permeated by a single set of parallel vertical fractures. Based on the anisotropic poroelasticity and perturbation theory, we first derive a linearized amplitude versus offset and azimuth approximation using the scattering function to decouple the fluid indicator and fracture parameters. Incorporating Bayes formula and convolution theory, we propose a feasible direct inversion approach in a Bayesian framework to obtain the direct estimations of model parameters, in which Cauchy and Gaussian distribution are used for the a priori information of model parameters and the likelihood function, respectively. We finally use the non‐linear iteratively reweighted least squares to solve the maximum a posteriori solutions of model parameters. The synthetic examples containing a moderate noise demonstrate the feasibility of the proposed approach, and the real data illustrates the stabilities of estimated fluid indicator and dry fracture parameters in gas‐saturated fractured porous rocks.
-
-
-
Amplitude variation with incident angle and azimuth inversion for Young's impedance, Poisson's ratio and fracture weaknesses in shale gas reservoirs
Authors Xinpeng Pan and Guangzhi ZhangABSTRACTBy analogy with P‐ and S‐wave impedances, the product of Young's modulus and density can be termed as Young's impedance, which indicates the rock lithology and brittleness of unconventional hydrocarbon reservoirs. Poisson's ratio is also an effective indicator of rock brittleness and fluid property of unconventional reservoirs, and fracture weaknesses indicate the fracture properties (fracturing intensity and fracture fillings) in fracture‐induced unconventional reservoirs. We aim to simultaneously estimate the Young's impedance, Poisson's ratio and fracture weaknesses from wide‐azimuth surface seismic data in a fracture‐induced shale gas reservoir, and use the horizontal transversely isotropic model to characterize the fractures. First, the linearized PP‐wave reflection coefficient in terms of Young's impedance, Poisson's ratio, density and fracture weaknesses is derived for the case of a weak‐contrast interface separating two weakly horizontal transversely isotropic media. In addition, an orthorhombic anisotropic case is also discussed in this paper. Then a Bayesian amplitude variation with incident angle and azimuth scheme with a model constraint is used to stably estimate Young's impedance, Poisson's ratio and fracture weaknesses with only PP‐wave azimuthal seismic data. The proposed approach is finally demonstrated on both synthetic and real data sets with reasonable results.
-
-
-
Rock physics modelling and inversion for saturation‐pressure changes in time‐lapse seismic studies
Authors Xiaozheng Lang and Dario GranaABSTRACTTime‐lapse seismic data are generally used to monitor the changes in dynamic reservoir properties such as fluid saturation and pore or effective pressure. Changes in saturation and pressure due to hydrocarbon production usually cause changes in the seismic velocities and as a consequence changes in seismic amplitudes and travel times. This work proposes a new rock physics model to describe the relation between saturation‐pressure changes and seismic changes and a probabilistic workflow to quantify the changes in saturation and pressure from time‐lapse seismic changes. In the first part of this work, we propose a new quadratic approximation of the rock physics model. The novelty of the proposed formulation is that the coefficients of the model parameters (i.e. the saturation‐pressure changes) are functions of the porosity, initial saturation and initial pressure. The improvements in the results of the forward model are shown through some illustrative examples. In the second part of the work, we present a Bayesian inversion approach for saturation‐pressure 4D inversion in which we adopt the new formulation of the rock physics approximation. The inversion results are validated using synthetic pseudo‐logs and a 3D reservoir model for CO2 sequestration.
-
-
-
Three‐dimensional electrical exploration methods for the mapping of polymetallic targets in Gansu Province, China
Authors Shengping Gong, Yabin Yang, Pinrong Lin, Wenli Wu, Caijun Zheng, Fusheng Shi, Xiaoping Wu, Aihua Weng, Guangzhi Zhang, Guanwen Gu and Yixin YeABSTRACTWe describe a three‐dimensional electrical exploration system, and its use in carrying out a three‐dimensional survey on a polymetallic deposit located in China's Gansu Province. A time‐domain‐induced polarization method, along with a controlled‐source audiomagnetotelluric method, was used to survey the same station points for comparison. Then, a three‐dimensional inversion based on the incomplete Gauss–Newton method (for time‐domain‐induced polarization) and a nonlinear conjugate gradient (for controlled‐source audiomagnetotelluric) method were developed and applied in order to model the data sets. Four subsurface targets were mapped by the three‐dimensional electrical exploration method. A strong correlation was found between the zones of the revealed ore‐bearing anomalies and the geological setting. Using the results of the three‐dimensional electrical exploration, we revealed the spatial distribution of the ore‐bearing stratum, as well as the relationship between the ore bodies and the stratum. The results potentially provided evidence for future interpretations of deep geological structures, along with evaluations of mineral resources.
-
-
-
A hybrid optimization scheme for self‐potential measurements due to multiple sheet‐like bodies in arbitrary 2D resistivity distributions
ABSTRACTSelf‐potential is a passive geophysical method that can be applied in a straightforward manner with minimum requirements in the field. Nonetheless, interpretation of self‐potential data is particularly challenging due to the inherited non‐uniqueness present in all potential methods. Incorporating information regarding the target of interest can facilitate interpretation and increase the reliability of the final output. In the current paper, a novel method for detecting multiple sheet‐like targets is presented. A numerical framework is initially described that simulates sheet‐like bodies in an arbitrary 2D resistivity distribution. A scattered field formulation based on finite differences is employed that allows the edges of the sheet to be independent of the grid geometry. A novel analytical solution for two‐layered models is derived and subsequently used to validate the accuracy of the proposed numerical scheme. Lastly, a hybrid optimization is proposed that couples linear least‐squares with particle‐swarm optimization in order to effectively locate the edges of multiple sheet‐like bodies. Through numerical and real data, it is proven that the hybrid optimization overcomes local minimal that occurs in complex resistivity distributions and converges substantially faster compared to traditional particle‐swarm optimization.
-
-
-
Improved imaging of a karst aquifer using focused source electromagnetic and differentially normalized method: a qualitative analysis
More LessABSTRACTGeophysical investigations using conventional techniques applied to groundwater exploration can often present strong limitations involving high financial costs, complex acquisition logistics and high ambiguity in results. Dispersion of the electric current flow, induced polarization) effects, cultural noises and shallow lateral heterogeneities represent the main problems faced by geoelectric methods in these types of surveys. Moreover, elements such as intrusions and mineralization at different depths may be responsible for signal attenuation as well as high resistivity in unsaturated zones and complex three‐dimensional formations or clayey zones cause variations in the electric current. The focused source electromagnetic and differentially normalized method approaches can help to solve some these issues. Aiming at a higher signal‐to‐noise ratio, the focused source electromagnetic method and approaches of the differentially normalized method, first applied to petroleum exploration, are tested on a groundwater target, in a karst environment sectioned by a diabase dyke. We performed the processing and analysis on real IP resistivity profiling data acquired with two‐way dipole‐dipole array, guided by magnetic data acquired on the same profile, mapping a diabase dyke. The inversion of focused source electromagnetic method/differentially normalized method was not performed, instead that we converted the induced polarization–resistivity data to a differential signal to qualitatively prove the presence of aquifer. Joint interpretation of focused source electromagnetic method curves and inverted two‐dimensional induced polarization–resistivity sections allowed for precise delineation of a conductive zone associated with the karst aquifer, le magnetics allowed for the definition of a neighbour dyke. The techniques have great potential in the aid of groundwater exploration, contributing substantially to the reduction of interpretation ambiguity. Focused source electromagnetic method/differentially normalized method/ approaches show that a simple linear combination of the conventional geoelectric data is able to remove the geological noise and provide the vertical focusing of the electric current.
-
Volumes & issues
-
Volume 72 (2023 - 2024)
-
Volume 71 (2022 - 2023)
-
Volume 70 (2021 - 2022)
-
Volume 69 (2021)
-
Volume 68 (2020)
-
Volume 67 (2019)
-
Volume 66 (2018)
-
Volume 65 (2017)
-
Volume 64 (2015 - 2016)
-
Volume 63 (2015)
-
Volume 62 (2014)
-
Volume 61 (2013)
-
Volume 60 (2012)
-
Volume 59 (2011)
-
Volume 58 (2010)
-
Volume 57 (2009)
-
Volume 56 (2008)
-
Volume 55 (2007)
-
Volume 54 (2006)
-
Volume 53 (2005)
-
Volume 52 (2004)
-
Volume 51 (2003)
-
Volume 50 (2002)
-
Volume 49 (2001)
-
Volume 48 (2000)
-
Volume 47 (1999)
-
Volume 46 (1998)
-
Volume 45 (1997)
-
Volume 44 (1996)
-
Volume 43 (1995)
-
Volume 42 (1994)
-
Volume 41 (1993)
-
Volume 40 (1992)
-
Volume 39 (1991)
-
Volume 38 (1990)
-
Volume 37 (1989)
-
Volume 36 (1988)
-
Volume 35 (1987)
-
Volume 34 (1986)
-
Volume 33 (1985)
-
Volume 32 (1984)
-
Volume 31 (1983)
-
Volume 30 (1982)
-
Volume 29 (1981)
-
Volume 28 (1980)
-
Volume 27 (1979)
-
Volume 26 (1978)
-
Volume 25 (1977)
-
Volume 24 (1976)
-
Volume 23 (1975)
-
Volume 22 (1974)
-
Volume 21 (1973)
-
Volume 20 (1972)
-
Volume 19 (1971)
-
Volume 18 (1970)
-
Volume 17 (1969)
-
Volume 16 (1968)
-
Volume 15 (1967)
-
Volume 14 (1966)
-
Volume 13 (1965)
-
Volume 12 (1964)
-
Volume 11 (1963)
-
Volume 10 (1962)
-
Volume 9 (1961)
-
Volume 8 (1960)
-
Volume 7 (1959)
-
Volume 6 (1958)
-
Volume 5 (1957)
-
Volume 4 (1956)
-
Volume 3 (1955)
-
Volume 2 (1954)
-
Volume 1 (1953)