- Home
- A-Z Publications
- Geophysical Prospecting
- Previous Issues
- Volume 68, Issue 4, 2020
Geophysical Prospecting - Volume 68, Issue 4, 2020
Volume 68, Issue 4, 2020
-
-
Amplitude‐compensated Laplacian filtering of reverse time migration and its application
Authors Renhu Yang, Xu Chang, Yun Ling, Ye E. Feng and Lijian QuABSTRACTReverse time migration is an advanced seismic migration imaging method. When the source wavefield and the receiver wavefield are cross‐correlated, the cross‐correlations of direct arrivals, backscattered waves and overturned waves will produce a lot of low‐frequency noise, which will mask the final imaging results. Laplacian filtering, as a common method to suppress low‐frequency noise, can adapt to any complex media, just adding a little computational cost. However, simple direct Laplacian filtering will destroy the characteristics of the useful signals. Therefore, the amplitude needs to be compensated before filtering when using the Laplacian filtering method. Zhang and Sun proposed an improved Laplacian filtering method and gave a simple calculation formula and explanation. This method can effectively suppress the low‐frequency noise in reverse time migration while retaining the useful signal characteristics, but lacks detailed and strict mathematical derivation. Therefore, this paper gives a detailed and rigorous mathematical derivation of the amplitude‐compensated Laplace filtering method from the point of view of amplitude‐preserved filtering. The source wavelet is used instead of the source wavefield to compensate amplitude, just adding a little calculation cost. Finally, the amplitude‐compensated Laplace filtering method is verified by two theoretical models and compared with the direct Laplacian filtering method.
-
-
-
Reverse time migration of ground‐penetrating radar with full wavefield decomposition based on the Hilbert transform
Authors Shichao Zhong, Yibo Wang, Yikang Zheng and Xu ChangABSTRACTThe conventional reverse time migration of ground‐penetrating radar data is implemented with the two‐way wave equation. The cross‐correlation result contains low‐frequency noise and false images caused by improper wave paths. To eliminate low‐frequency noise and improve the quality of the migration image, we propose to separate the left‐up‐going, left‐down‐going, right‐up‐going and right‐down‐going wavefield components in the forward‐ and backward‐propagated wavefields based on the Hilbert transform. By applying the reverse time migration of ground‐penetrating radar data with full wavefield decomposition based on the Hilbert transform, we obtain the reverse time migration images of different wavefield components and combine correct imaging conditions to generate complete migration images. The proposed method is tested on the synthetic ground‐penetrating radar data of a tilt‐interface model and a complex model. The migration results show that the imaging condition of different wavefield components can highlight the desired structures. We further discuss the reasons for incomplete images by reverse time migration with partial wavefields. Compared with the conventional reverse time migration methods for ground‐penetrating radar data, low‐frequency noise can be eliminated in images generated by the reverse time migration method with full wavefield decomposition based on the Hilbert transform.
-
-
-
Coherence algorithm with a high‐resolution time–time transform and feature matrix for seismic data
Authors Fengyuan Sun, Jinghuai Gao, Bing Zhang and Naihao LiuABSTRACTTraditional coherence algorithms are often based on the assumption that seismic traces are stationary and Gaussian. However, seismic traces are actually non‐stationary and non‐Gaussian. A constant time window and the canonical correlation analysis in traditional coherence algorithms are not optimal for non‐stationary seismic traces and cannot describe the similarity between adjacent seismic traces in detail. To overcome this problem, a new coherence algorithm using the high‐resolution time–time transform and the feature matrix is designed. The high‐resolution time–time transform used to replace the constant time window can produce a frequency‐dependent time local series to analyse non‐stationary seismic traces. The feature matrix, constructed by the frequency‐dependent time local series and the related local gradients, defines a new correlation metric that enhances more details of the geological discontinuities in seismic images than does the canonical correlation analysis. Additionally, the Riemannian metric is introduced for related calculations because the feature matrices are not defined in a Euclidean space but rather in a manifold space. Application to field data illustrates that the proposed method reveals more details of structural and stratigraphic features.
-
-
-
Triplications for the converted wave in transversely isotropic media with a tilted symmetry axis
Authors Shibo Xu, Alexey Stovas and Hitoshi MikadaABSTRACTIn seismic data processing, serious problems could be caused by the existence of triplication and need to be treated properly for tomography and other inversion methods. The triplication in transversely isotropic medium with a vertical symmetry axis has been well studied and concluded that the triplicated traveltime only occurs for S wave and there is no triplication for P and converted PS waves since the P wave convexity slowness always compensates the S wave slowness concavity. Compared with the vertical symmetry axis model, the research of the triplication in transversely isotropic medium with a tilted symmetry axis is still keeping blank. In order to analyse the triplication for the converted wave in the tilted symmetry axis model, we examine the traveltime of the triplication from the curvature of averaged P and S wave slowness. Three models are defined and tested in the numerical examples to illustrate the behaviour of the tilted symmetry axis model for the triplicated traveltime with the change of the rotation angle. Since the orientation of an interface is related to the orientation of the symmetry axis, the triplicated traveltime is encountered for the converted wave in the tilted symmetry axis model assuming interfaces to be planar and horizontal. The triplicated region is influenced by the place and level of the concave curvature of the P and S wave slowness.
-
-
-
Results of a walk‐above vertical seismic profiling survey acquired at the Thônex‐01 geothermal well (Switzerland) to delineate fractured carbonate formations for geothermal development
ABSTRACTThis paper deals with the design, acquisition and processing of vertical seismic profiling data collected in 2016 at the directional Thônex‐01 geothermal well, drilled in 1993 in the Geneva Canton (Switzerland). The aim of the study was to obtain additional information from existing‐abandoned well for downhole geophysical characterization of potential geothermal reservoirs in the Geneva area. The main results obtained by this study allowed an improved subsurface image by multi‐offset acquisitions near and around the deviated well. Specifically we were able to (1) improve the velocity model of the Geneva Basin, which, at present, is only constrained by one unique deep exploration well, located at the south‐west of the basin, far away from the study area; (2) image the transition between the Tertiary Molasse siliciclastic sediments and the Mesozoic carbonates; (3) define an acquisition strategy for future vertical seismic profiling surveys in support of the geothermal exploration campaigns in the area. This study has demonstrated the value of vertical seismic profiling method for the refinement of the subsurface geology at local scale, and its usefulness to assist the planning of second well associated with geothermal doublets.
-
-
-
The effectiveness of a pseudo‐inverse extended Born operator to handle lateral heterogeneity for imaging and velocity analysis applications
Authors Abdullah Alali, Bingbing Sun and Tariq AlkhalifahABSTRACTWave equation–based migration velocity analysis techniques aim to construct a kinematically accurate velocity model for imaging or as an initial model for full waveform inversion applications. The most popular wave equation–based migration velocity analysis method is differential semblance optimization, where the velocity model is iteratively updated by minimizing the unfocused energy in an extended image volume. However, differential semblance optimization suffers from artefacts, courtesy of the adjoint operator used in imaging, leading to poor convergence. Recent findings show that true amplitude imaging plays a significant role in enhancing the differential semblance optimization's gradient and reducing the artefacts. Here, we focus on a pseudo‐inverse operator to the horizontally extended Born as a true amplitude imaging operator. For laterally inhomogeneous models, the operator required a derivative with respect to a vertical shift. Extending the image vertically to evaluate such a derivative is costly and impractical. The inverse operator can be simplified in laterally homogeneous models. We derive an extension of the approach to apply the full inverse formula and evaluate the derivative efficiently. We simplified the implementation by applying the derivative to the imaging condition and utilize the relationship between the source and receiver wavefields and the vertical shift. Specifically, we verify the effectiveness of the approach using the Marmousi model and show that the term required for the lateral inhomogeneity treatment has a relatively small impact on the results for many cases. We then apply the operator in differential semblance optimization and invert for an accurate macro‐velocity model, which can serve as an initial velocity model for full waveform inversion.
-
-
-
Simultaneous joint migration inversion for high‐resolution imaging/inversion of time‐lapse seismic datasets
Authors Shan Qu and Dirk Jacob VerschuurABSTRACTThe current time‐lapse practice is to exactly repeat well‐sampled acquisition geometries to mitigate acquisition effects on the time‐lapse differences. In order to relax the rigid requirements on acquisition effects, we propose simultaneous joint migration inversion as an effective time‐lapse tool for reservoir monitoring, which combines a joint time‐lapse data processing strategy with the joint migration inversion method. Joint migration inversion is a full‐wavefield inversion method that explains the measured reflection data using a parameterization in terms of reflectivity and propagation velocity. Both the inversion process inside the imaging/inversion scheme and the extra illumination provided by including multiples in joint migration inversion makes the obtained velocity and reflectivity operator largely independent of the utilized acquisition geometry and, thereby, relaxes the strong requirement of non‐repeatability during the monitoring. Because simultaneous joint migration inversion inverts for all datasets simultaneously and utilizes various constraints on the estimated reflectivities and velocity, the obtained time‐lapse differences have much higher accuracy compared to inverting each dataset separately. It allows the baseline and monitor parameters to communicate with each other dynamically during inversion via a user‐defined spatial weighting operator. In order to get more localized time‐lapse velocity differences, we further extend the regular simultaneous joint migration inversion to a robust high‐resolution simultaneous joint migration inversion process using the time‐lapse reflectivity difference as an extra constraint for the velocity estimation during inversion. This constraint makes a link between the reflectivity‐ and the velocity difference by exploiting the relationship between them. We demonstrate the feasibility of the proposed method with a highly realistic synthetic model based on the Grane field offshore Norway and a time‐lapse field dataset from the Troll Field.
-
-
-
An efficiency‐improved genetic algorithm and its application on multimodal functions and a 2D common reflection surface stacking problem
Authors Yenni Paloma Villa Acuna and Yimin SunABSTRACTAlthough Genetic Algorithms have found many successful applications in the field of exploration geophysics, the convergence speed remains a big challenge as Genetic Algorithms usually require a huge amount of fitness function evaluations. In this paper, we propose an efficiency‐improved Genetic Algorithm, which has both a good global search capability and a good local search capability, and is also capable of robustly handling the premature convergence challenge commonly seen in linear and directed non‐linear optimization methods. In our new genetic algorithm, the global search capability is performed via a modified island model, while the local search capability is provided by a novel self‐adaptive differential evolution fine tuning scheme. Premature convergence is dealt with via a local exhaustive search method. We first demonstrate the much improved convergence speed of this efficiency‐improved Genetic Algorithm over that of our previously proposed advanced Genetic Algorithm on several multimodal functions. We further demonstrate the effectiveness of our efficiency‐improved Genetic Algorithm by applying it to a two‐dimensional common reflection surface stacking problem, which is a highly nonlinear geophysical optimization problem, to obtain very encouraging results.
-
-
-
Scattering pattern analysis and true‐amplitude generalized Radon transform migration for acoustic transversely isotropic media with a vertical axis of symmetry
Authors Quan Liang, Wei Ouyang, Weijian Mao and Shijun ChengABSTRACTIn multi‐parameter ray‐based anisotropic migration/inversion, it is essential that we have an understanding of the scattering mechanism corresponding to parameter perturbations. Because the complex nonlinearity in the anisotropic inversion problem is intractable, the construction of true‐amplitude linearized migration/inversion procedures is needed and important. By using the acoustic medium assumption for transversely isotropic media with a vertical axis of symmetry and representing the anisotropy with P‐wave normal moveout velocity, Thomsen parameter δ and anelliptic parameter η, we formalize the linearized inverse scattering problem for three‐dimensional pseudo‐acoustic equations. Deploying the single‐scattering approximation and an elliptically anisotropic background introduces a new linear integral operator that connects the discontinuous perturbation parameters with the multi‐shot/multi‐offset P‐wave scattered data. We further apply the high‐frequency asymptotic Green's function and its derivatives to the integral operator, and then the scattering pattern of each perturbation parameter can be explicitly presented. By naturally establishing a connection to generalized Radon transform, the pseudo‐inverse of the integral operator can be solved by the generalized Radon transform inversion. In consideration of the structure of this pseudo‐inverse operator, the computational implementation is done pointwise by shooting a fan of rays from the target imaging area towards the acquisition system. Results from two‐dimensional numerical tests show amplitude‐preserving images with high quality.
-
-
-
A high accurate automated first‐break picking method for seismic records from high‐density acquisition in areas with a complex surface
Authors Yinpo Xu, Cheng Yin, Xuefeng Zou, Yudong Ni, Yingjie Pan and Liangjun XuABSTRACTAs the application of high‐density high‐efficiency acquisition technology becomes more and more wide, the areas with complex surface conditions gradually become target exploration areas, and the first‐break picking work of massive low signal‐to‐noise ratio data is a big challenge. The traditional method spends a lot of manpower and time to interactively pick first breaks, a large amount of interactive work affects the accuracy and efficiency of picking. In order to overcome the shortcoming that traditional methods have weak anti‐noise to low signal‐to‐noise ratio primary wave, this paper proposes a high accurate automated first‐break picking method for low signal‐to‐noise ratio primary wave from high‐density acquisition in areas with a complex surface. Firstly, this method determines first‐break time window using multi‐azimuth spatial interpolation technology; then it uses the improved clustering algorithm to initially pick first breaks and then perform multi‐angle comprehensive quality evaluation to first breaks according to the following sequence: ‘single trace → spread → single shot → multiple shots’ to identify the abnormal first breaks; finally it determines the optimal path through the constructed evaluation function and using the ant colony algorithm to correct abnormal first breaks. Multi‐azimuth time window spatial interpolation technology provides the base for accurately picking first‐break time; the clustering algorithm can effectively improve the picking accuracy rate of low signal‐to‐noise ratio primary waves; the multi‐angle comprehensive quality evaluation can accurately and effectively eliminate abnormal first breaks; the ant colony algorithm can effectively improve the correction quality of low signal‐to‐noise ratio abnormal first breaks. By example analysis and comparing with the commonly used Akaike Information Criterion method, the automated first‐break picking theory and technology studied in this paper has high picking accuracy and the ability to stably process low signal‐to‐noise ratio seismic data, has a significant effect on seismic records from high‐density acquisition in areas with a complex surface and can meet the requirements of accuracy and efficiency for massive data near‐surface modelling and statics calculation.
-
-
-
A nearly analytic symplectic partitioned Runge–Kutta method based on a locally one‐dimensional technique for solving two‐dimensional acoustic wave equations
Authors Chol Sim, Chunyou Sun and Nam YunABSTRACTIn this paper, we develop a new nearly analytic symplectic partitioned Runge–Kutta method based on locally one‐dimensional technique for numerically solving two‐dimensional acoustic wave equations. We first split two‐dimensional acoustic wave equation into the local one‐dimensional equations and transform each of the split equations into a Hamiltonian system. Then, we use both a nearly analytic discrete operator and a central difference operator to approximate the high‐order spatial differential operators, which implies the symmetry of the discretized spatial differential operators, and we employ the partitioned second‐order symplectic Runge–Kutta method to numerically solve the resulted semi‐discrete Hamiltonian ordinary differential equations, which results in fully discretized scheme is symplectic unlike conventional nearly analytic symplectic partitioned Runge–Kutta methods. Theoretical analyses show that the nearly analytic symplectic partitioned Runge–Kutta method based on locally one‐dimensional technique exhibits great higher stability limits and less numerical dispersion than the nearly analytic symplectic partitioned Runge–Kutta method. Numerical experiments are conducted to verify advantages of the nearly analytic symplectic partitioned Runge–Kutta method based on locally one‐dimensional technique, such as their computational efficiency, stability, numerical dispersion and long‐term calculation capability.
-
-
-
Separation of multi‐mode surface waves by supervised machine learning methods
Authors Jing Li, Yuqing Chen and Gerard T. SchusterABSTRACTLogistic regression, neural networks and support vector machines are tested for their effectiveness in isolating surface waves in seismic shot records. To distinguish surface waves from other arrivals, we train the algorithms on three distinguishing features of surface‐wave dispersion curves in the domain: spectrum coherency of the trace's magnitude spectrum, local dip and the frequency range for a fixed wavenumber k in the spectrum. Numerical tests on synthetic data show that the kernel‐based support vector machines algorithm gives the highest accuracy in predicting the surface‐wave window in the domain compared to neural networks and logistic regression. This window is also used to automatically pick the fundamental dispersion curve. The other two methods correctly pick the low‐frequency part of the dispersion curve but fail at higher frequencies where there is interference with higher‐order modes.
-
-
-
Green's function for three‐dimensional elastic wave equation with a moving point source on the free surface with applications
Authors Jing‐Bo Chen and Jian CaoABSTRACTHigh‐speed train seismology has come into being recently. This new kind of seismology uses a high‐speed train as a repeatable moving seismic source. Therefore, Green's function for a moving source is needed to make theoretical studies of the high‐speed train seismology. Green's function for three‐dimensional elastic wave equation with a moving point source on the free surface is derived. It involves a line integral of the Green's function for a fixed point source with different positions and corresponding time delays. We give a rigorous mathematical proof of this Green's function. According to the principle of linear superposition, we have also obtained the Green's function for a group of moving sources which can be regarded as a model of a traveling high‐speed train. Based on a temporal convolution, an analytical formula for other moving sources is also given. In terms of a moving Gaussian source, we deal with the issue of numerical calculations of the analytical formula. Applications to modelling of a traveling high‐speed train are presented. We have considered both the land case and the bridge case for a traveling high‐speed train. The theoretical seismograms show different waveform features for these two cases.
-
-
-
A hybrid approach to compute seismic travel times in three‐dimensional tetrahedral meshes
Authors Maher Nasr, Bernard Giroux and J. Christian DupuisABSTRACTWe propose an optimized method to compute travel times for seismic inversion problems. It is a hybrid method combining several approaches to deal with travel time computation accuracy in unstructured meshes based on tetrahedral elementary cells. As in the linear travel time interpolation method, the proposed approach computes travel times using seismic ray paths. The method operates in two sequential steps: At a first stage, travel times are computed for all nodes of the mesh using a modified version of the shortest path method. The difference with the standard version is that additional secondary nodes (called tertiary nodes) are added temporarily around seismic sources in order to improve accuracy with a reasonable increase in computational cost. During the second step, the steepest travel time gradient method is used to trace back ray paths for each source–receiver pair. Travel times at each receiver are then recomputed using slowness values at the intersection points between the ray path and the traversed cells. A number of numerical tests with an array of different velocity models, mesh resolutions and mesh topologies have been carried out. These tests showed that an average relative error in the order of 0.1% can be achieved at a computational cost that is suitable for travel time inversion.
-
-
-
Establishing an effective offset‐dependent velocity model by ray tracing and its application in Kirchhoff time migration
Authors Guiting Chen, Zhenli Wang, Oluwatosin John Rotimi and Suyang ZhangABSTRACTConventional Kirchhoff prestack time migration based on the hyperbolic moveout can cause ambiguity in laterally inhomogeneous media, because the root mean square velocity corresponds to a one‐dimensional model under the horizontal layer assumption; it does not include the lateral variations. The shot/receiver configuration with different offsets and azimuths should adopt different migration velocities as they contribute to a single image point. Therefore, we propose to use an offset‐vector to describe the lateral variations through an offset‐dependent velocity corresponding to the difference in offset from surface points to the image point. The offset‐vector is decomposed into orthogonal directions along the in‐line and cross‐line directions so that the single velocity can be expressed as a series of actual velocities. We use a simple Snell's law‐based ray tracing to calculate the travel time recorded at the image point and convert the travel time to an equivalent velocity corresponding to a pseudo‐straight ray. The double‐square‐root equation using such an equivalent velocity in the offset‐vector domain is non‐hyperbolic and asymmetrical, which improves the accuracy of the migration. Numerical examples using the Marmousi model and a wide azimuth field data show that the proposed method can achieve reasonable accuracy and significantly enhances the imaging of complex structures.
-
-
-
Estimation of the seismic wavelet through homomorphic deconvolution and well log data: application on well‐to‐seismic tie procedure
ABSTRACTWavelet estimation and well‐tie procedures are important tasks in seismic processing and interpretation. Deconvolutional statistical methods to estimate the proper wavelet, in general, are based on the assumptions of the classical convolutional model, which implies a random process reflectivity and a minimum‐phase wavelet. The homomorphic deconvolution, however, does not take these premises into account. In this work, we propose an approach to estimate the seismic wavelet using the advantages of the homomorphic deconvolution and the deterministic estimation of the wavelet, which uses both seismic and well log data. The feasibility of this approach is verified on well‐to‐seismic tie from a real data set from Viking Graben Field, North Sea, Norway. The results show that the wavelet estimated through this methodology produced a higher quality well tie when compared to methods of estimation of the wavelet that consider the classical assumptions of the convolutional model.
-
-
-
Elastic full waveform inversion with probabilistic petrophysical clustering
Authors Odette Aragao and Paul SavaABSTRACTFull waveform inversion aims to use all information provided by seismic data to deliver high‐resolution models of subsurface parameters. However, multiparameter full waveform inversion suffers from an inherent trade‐off between parameters and from ill‐posedness due to the highly non‐linear nature of full waveform inversion. Also, the models recovered using elastic full waveform inversion are subject to local minima if the initial models are far from the optimal solution. In addition, an objective function purely based on the misfit between recorded and modelled data may honour the seismic data, but disregard the geological context. Hence, the inverted models may be geologically inconsistent, and not represent feasible lithological units. We propose that all the aforementioned difficulties can be alleviated by explicitly incorporating petrophysical information into the inversion through a penalty function based on multiple probability density functions, where each probability density function represents a different lithology with distinct properties. We treat lithological units as clusters and use unsupervised K‐means clustering to separate the petrophysical information into different units of distinct lithologies that are not easily distinguishable. Through several synthetic examples, we demonstrate that the proposed framework leads full waveform inversion to elastic models that are superior to models obtained either without incorporating petrophysical information, or with a probabilistic penalty function based on a single probability density function.
-
-
-
Well data mapping validation
Authors D. Krilov and M. KrylovaABSTRACTSome key geologic data can only be obtained from wells. This article provides original technology concerning the geological mapping validation of complex anomalous zones based on information obtained from irregular or sparse drilling. Only complex non‐unique cases of geologic characteristics delineation (i.e. characterized by multiple solutions) are considered here. Representative groups of multiple mapping solutions (i.e. essential variations in anomaly zones delineation in terms of size and location) are analysed in order to optimize and verify the mapping parameters. Minor interpolation discrepancies not influencing the geologic interpretation, for example, resulting from specificities of different interpolation techniques or insignificant mapping procedure corrections, are neglected. The mapping error can be quantitatively assessed by the cross validation error analysis. Maps can be optimized by the proper choice of the interpolation algorithm and data filtering. Using trial‐and‐error schema, interpolation and filtering options are tested and mapping errors are compared. Single assessments or selective tests of the mapping accuracy often cannot lead to the best solution. The optimization procedure can be used for mapping verification, accuracy assessment and mapping parameter optimization. It may also be used to find an optimal balance between the mapping accuracy and the resolution by interpolants (i.e. initial well data) test processing. The interpolation algorithms used here are most technologically fit for a check comparison – not as a substitution of the existing mapping technologies. After positive mapping validation, a common procedure of interpolation can be started. Franke's theoretical work on scattered data interpolation and Krilov's smart averaging interpolation tests were used here to substantiate the well data mapping validation technology. The authors tried not to delve deep into the theory and interpolation testing details using due references and concentrating on the technology procedure and algorithm description.
-
-
-
Incorporating known petrophysical model in the seismic full‐waveform inversion using the Gramian constraint
Authors M. Malovichko, N. Khokhlov, N. Yavich and M. S. ZhdanovABSTRACTIn this paper, we develop a general approach to integrating petrophysical models in three‐dimensional seismic full‐waveform inversion based on the Gramian constraints. In the framework of this approach, we present an example of the frequency‐domain P‐wave velocity inversion guided by an electrical conductivity model. In order to introduce a coupling between the two models, we minimize the corresponding Gramian functional, which is included in the Tikhonov parametric functional. We demonstrate that in the case of a single‐physics inversion guided by a model of different physical type, the general expressions of the Gramian functional and its gradients become simple and easy to program. We also prove that the Gramian functional has a non‐negative quadratic form, so it can be easily incorporated in a standard gradient‐based minimization scheme. The developed new approach of seismic inversion guided by the known petrophysical model has been validated by three‐dimensional inversion of synthetic seismic data generated for a realistic three‐dimensional model of the subsurface.
-
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)