- Home
- A-Z Publications
- Geophysical Prospecting
- Previous Issues
- Volume 69, Issue 1, 2021
Geophysical Prospecting - Volume 69, Issue 1, 2021
Volume 69, Issue 1, 2021
-
-
Eigenrays in 3D heterogeneous anisotropic media, Part I: Kinematics
Authors Zvi Koren and Igor RavveABSTRACTWe present a new ray bending approach, referred to as the Eigenray method, for solving two‐point boundary‐value kinematic and dynamic ray tracing problems in 3D smooth heterogeneous general anisotropic elastic media. The proposed Eigenray method is aimed to provide reliable stationary ray path solutions and their dynamic characteristics, in cases where conventional initial‐value ray shooting methods, followed by numerical convergence techniques, become challenging. The kinematic ray bending solution corresponds to the vanishing first traveltime variation, leading to a stationary path between two fixed endpoints (Fermat's principle), and is governed by the nonlinear second‐order Euler–Lagrange equation. The solution is based on a finite‐element approach, applying the weak formulation that reduces the Euler–Lagrange second‐order ordinary differential equation to the first‐order weighted‐residual nonlinear algebraic equation set. For the kinematic finite‐element problem, the degrees of freedom are discretized nodal locations and directions along the ray trajectory, where the values between the nodes are accurately and naturally defined with the Hermite polynomial interpolation. The target function to be minimized includes two essential penalty (constraint) terms, related to the distribution of the nodes along the path and to the normalization of the ray direction. We distinguish between two target functions triggered by the two possible types of stationary rays: a minimum traveltime and a saddle‐point solution (due to caustics). The minimization process involves the computation of the global (all‐node) traveltime gradient vector and the traveltime Hessian matrix. The traveltime Hessian is used for the minimization process, analysing the type of the stationary ray, and for computing the geometric spreading of the entire resolved stationary ray path. The latter, however, is not a replacement for the dynamic ray tracing solution, since it does not deliver the geometric spreading for intermediate points along the ray, nor the analysis of caustics. Finally, we demonstrate the efficiency and accuracy of the proposed method along three canonical examples.
-
-
-
Eigenrays in 3D heterogeneous anisotropic media, Part II: Dynamics
Authors Igor Ravve and Zvi KorenABSTRACTThis paper is the second in a sequel of two papers and dedicated to the computation of paraxial rays and dynamic characteristics along the stationary rays obtained in the first paper. We start by formulating the linear, second‐order, Jacobi dynamic ray tracing equation. We then apply a similar finite‐element solver, as used for the kinematic ray tracing, to compute the dynamic characteristics between the source and any point along the ray. The dynamic characteristics in our study include the relative geometric spreading and the phase correction due to caustics (i.e. the amplitude and the phase of the asymptotic form of the Green's function for waves propagating in 3D heterogeneous general anisotropic elastic media). The basic solution of the Jacobi equation is a shift vector of a paraxial ray in the plane normal to the ray direction at each point along the central ray. A general paraxial ray is defined by a linear combination of up to four basic vector solutions, each corresponds to specific initial conditions related to the ray coordinates at the source. We define the four basic solutions with two pairs of initial condition sets: point–source and plane‐wave. For the proposed point–source ray coordinates and initial conditions, we derive the ray Jacobian and relate it to the relative geometric spreading for general anisotropy. Finally, we introduce a new dynamic parameter, similar to the endpoint complexity factor, presented in the first paper, used to define the measure of complexity of the propagated wave/ray phenomena. The new weighted propagation complexity accounts for the normalized relative geometric spreading not only at the receiver point, but along the whole stationary ray path. We propose a criterion based on this parameter as a qualifying factor associated with the given ray solution. To demonstrate the implementation of the proposed method, we use several isotropic and anisotropic benchmark models. For all the examples, we first compute the stationary ray paths, and then compute the geometric spreading and analyse these trajectories for possible caustics. Our primary aim is to emphasize the advantages, transparency and simplicity of the proposed approach.
-
-
-
Geometrical characteristics of phase and group velocity surfaces in anisotropic media
Authors Alexey Stovas, Yuriy Roganov and Vyacheslav RoganovABSTRACTThe phase and group velocity surfaces are essential for wave propagation in anisotropic media. These surfaces have certain features that, especially, for shear waves result in complications for modelling and inversion of recorded wavefields. To analyse wave propagation in an anisotropic model, it is important to identify these features in both the phase and group domains. We propose few characteristics for this analysis: the energy flux angle, decomposed in the polar and azimuth angle correction angles and enhancement factor, which is able to characterize both singularity points and triplication zones. The very simple equation that controls the triplications is derived in the phase domain. The proposed characteristics are illustrated for elastic and acoustic anisotropic models of different symmetry classes.
-
-
-
Time evolution of the first‐order linear acoustic/elastic wave equation using Lie product formula and Taylor expansion
Authors Edvaldo S. Araujo and Reynam C. PestanaABSTRACTWe propose a new numerical solution to the first‐order linear acoustic/elastic wave equation. This numerical solution is based on the analytic solution of the linear acoustic/elastic wave equation and uses the Lie product formula, where the time evolution operator of the analytic solution is written as a product of exponential matrices where each exponential matrix term is then approximated by Taylor series expansion. Initially, we check the proposed approach numerically and then demonstrate that it is more accurate to apply a Taylor expansion for the exponential function identity rather than the exponential function itself. The numerical solution formulated employs a recursive procedure and also incorporates the split perfectly matched layer boundary condition. Thus, our scheme can be used to extrapolate wavefields in a stable manner with even larger time‐steps than traditional finite‐difference schemes. This new numerical solution is examined through the comparison of the solution of full acoustic wave equation using the Chebyshev expansion approach for the matrix exponential term. Moreover, to demonstrate the efficiency and applicability of our proposed solution, seismic modelling results of three geological models are presented and the processing time for each model is compared with the computing time taking by the Chebyshev expansion method. We also present the result of seismic modelling using the scheme based in Lie product formula and Taylor series expansion for the first‐order linear elastic wave equation in vertical transversely isotropic and tilted transversely isotropic media as well. Finally, a post‐stack migration results are also shown using the proposed method.
-
-
-
Two‐dimensional dictionary learning for suppressing random seismic noise
Authors Kunhong Li, Zhining Liu, Bin She, Hanpeng Cai, Yaojun Wang and Guangmin HuABSTRACTDictionary learning is a successful method for random seismic noise attenuation that has been proven by some scholars. Dictionary learning–based techniques aim to learn a set of common bases called dictionaries from given noised seismic data. Then, the denoising process will be performed by assuming a sparse representation on each small local patch of the seismic data over the learned dictionary. The local patches that are extracted from the seismic section are essentially two‐dimensional matrices. However, for the sake of simplicity, almost all of the existing dictionary learning methods just convert each two‐dimensional patch into a one‐dimensional vector. In doing this, the geometric structure information of the raw data will be revealed, leading to low capability in the reconstruction of seismic structures, such as faults and dip events. In this paper, we propose a two‐dimensional dictionary learning method for the seismic denoising problem. Unlike other dictionary learning–based methods, the proposed method represents the two‐dimensional patches directly to avoid the conversion process, and thus reserves the important structure information for a better reconstruction. Our method first learns a two‐dimensional dictionary from the noisy seismic patches. Then, we use the two‐dimensional dictionary to sparsely represent all of the noisy two‐dimensional patches to obtain clean patches. Finally, the clean patches are patched back to generate a denoised seismic section. The proposed method is compared with the other three denoising methods, including FX‐decon, curvelet and one‐dimensional learning method. The results demonstrate that our method has better denoising performance in terms of signal‐to‐noise ratio, fault and amplitude preservation.
-
-
-
Erratic noise suppression using iterative structure‐oriented space‐varying median filtering with sparsity constraint
Authors Guangtan Huang, Min Bai, Qiang Zhao, Wei Chen and Yangkang ChenABSTRACTErratic noise often has high amplitudes and a non‐Gaussian distribution. Least‐squares–based approaches therefore are not optimal. This can be handled better with non–least‐squares approaches, for example based on Huber norm which is computationally expensive. An alternative method has been published which involves transforming the data with erratic noise to pseudodata that have Gaussian distributed noise. It can then be attenuated using traditional least‐squares approaches. This alternative method has previously been used in combination with a curvelet transform in an iterative scheme. In this paper, we introduce a median‐filtering step in this iterative scheme. The median filter is applied following the slope direction of the seismic data to maximally preserve the energy of useful signals. The new method can suppress stronger erratic noise compared with the previous iterative method, and can better deal with random noise compared with the single‐step implementation of the median filter. We apply the proposed robust denoising algorithm to a synthetic dataset and two field data examples and demonstrate its advantages over three different noise attenuation algorithms.
-
-
-
Passive seismic data primary estimation and noise removal via focal‐denoising closed‐loop surface‐related multiple elimination based on 3D L1‐norm sparse inversion
Authors Tiexing Wang, Deli Wang, Jing Sun and Bin HuABSTRACTPassive seismic has recently attracted a great deal of attention because non‐artificial source is used in subsurface imaging. The utilization of passive source is low cost compared with artificial‐source exploration. In general, constructing virtual shot gathers by using cross‐correlation is a preliminary step in passive seismic data processing, which provides the basis for applying conventional seismic processing methods. However, the subsurface structure is not uniformly illuminated by passive sources, which leads to that the ray path of passive seismic does not fit the hyperbolic hypothesis. Thereby, travel time is incorrect in the virtual shot gathers. Besides, the cross‐correlation results are contaminated by incoherent noise since the passive sources are always natural. Such noise is kinematically similar to seismic events and challenging to be attenuated, which will inevitably reduce the accuracy in the subsequent process. Although primary estimation for transient‐source seismic data has already been proposed, it is not feasible to noise‐source seismic data due to the incoherent noise. To overcome the above problems, we proposed to combine focal transform and local similarity into a highly integrated operator and then added it into the closed‐loop surface‐related multiple elimination based on the 3D L1‐norm sparse inversion framework. Results proved that the method was capable of reliably estimating noise‐free primaries and correcting travel time at far offsets for a foresaid virtual shot gathers in a simultaneous closed‐loop inversion manner.
-
-
-
Calculating complex‐valued P‐wave first‐arrival traveltimes in attenuative vertical transversely isotropic media with an irregular surface
Authors Mengxiu Wang, Jingyi Chen, Zhencong Zhao and Wenyi HuABSTRACTThe complex‐valued first‐arrival traveltime can be used to describe the properties of both velocity and attenuation as seismic waves propagate in attenuative elastic media. The real part of the complex‐valued traveltime corresponds to phase arrival and the imaginary part is associated with the amplitude decay due to energy absorption. The eikonal equation for attenuative vertical transversely isotropic media discretized with rectangular grids has been proven effective and precise to calculate the complex‐valued traveltime, but less accurate and efficient for irregular models. By using the perturbation method, the complex‐valued eikonal equation can be decomposed into two real‐valued equations, namely the zeroth‐ and first‐order traveltime governing equations. Here, we first present the topography‐dependent zeroth‐ and first‐order governing equations for attenuative VTI media, which are obtained by using the coordinate transformation from the Cartesian coordinates to the curvilinear coordinates. Then, we apply the Lax–Friedrichs sweeping method for solving the topography‐dependent traveltime governing equations in order to approximate the viscosity solutions, namely the real and imaginary parts of the complex‐valued traveltime. Several numerical tests demonstrate that the proposed scheme is efficient and accurate in calculating the complex‐valued P‐wave first‐arrival traveltime in attenuative VTI media with an irregular surface.
-
-
-
Optimized experimental design for seismic full waveform inversion: A computationally efficient method including a flexible implementation of acquisition costs
Authors Valérie Krampe, Pascal Edme and Hansruedi MaurerABSTRACTOptimized experimental design aims at reducing the cost of a seismic survey by identifying the optimal locations and amounts of sources and receivers. While the acquisition design in the context of seismic imaging applies criteria like fold, offset and spatial sampling, different attributes such as the sensitivity kernels are more relevant for seismic full waveform inversion. An ideal measure to quantify the goodness of an acquisition design relies on the eigenvalue spectrum of the approximate Hessian matrix, but this technique is computationally too expensive for practical use. A more affordable goodness measure has been proposed in the past, but we demonstrate that this measure is inappropriate for target‐oriented optimized experimental design. To address those issues, we derived a sequential receiver‐based procedure using a goodness measure based on the determinant of the approximate Hessian matrix. We show with numerical tests that it efficiently provides an optimized design for target‐oriented as well as for extensive full waveform inversion. This design allows a better reconstruction of the subsurface than an evenly spaced acquisition geometry. Furthermore, the optimization algorithm itself can easily be parallelized, therefore making it attractive for applications to large‐scale three‐dimensional surveys. In addition, our algorithm is able to incorporate variable costs, representing any kind of acquisition‐related costs, for every individual source location. The combined optimization with respect to the information content of sources and to the true cost will allow a more comprehensive and realistic survey planning and has a high potential for further applications.
-
-
-
Mapping of underground cavities by the passive seismic standing waves method: the case study of Barsukovskaya cave (Novosibirsk region, Russia)
Authors Konstantin V. Fedin, Yury I. Kolesnikov and Luckymore NgomayezweABSTRACTThis article presents the results of mapping a karst cave by the passive seismic standing waves method. Barsukovskaya cave is located about 100 km southeast of the city of Novosibirsk (Russia). The total length of the cave's passages and grottoes is estimated at about 200 m, the maximum depth from the earth's surface is about 19 m. The method for studying underground cavities used is based on the effect of the generation of standing waves by microtremor in the space between the earth's surface and the cave roof. The accumulation of amplitude spectra of a large number of microtremor records makes it possible to determine the frequencies of the first few modes of these waves. Areal passive seismic survey on the earth's surface above the cave made it possible to construct a map of the lowest mode frequency distribution over the cave roof. Since no standing waves were observed at other points, this map reflects the cave structure in plan, which confirms the comparison with the cave diagram drawn up earlier by one of the speleologists. A schematic map of the depth of the cave roof was constructed using the longitudinal wave velocity Vp = 3120 m/s determined by the rock samples selected near the entrance to the cave. This map at a qualitative level also agrees with the data of speleologists, which indicate that the cave, on average, gradually becomes deeper from the entrance to its dead‐end branches. The shallower depths in comparison with the data of speleologists are apparently explained by a very low estimate of the velocity determined from a rock sample taken near the entrance to the cave. The reliability of the obtained cave mapping results is confirmed by the numerical simulation results using the finite‐element method.
-
-
-
Anisotropic attenuation of stratified viscoelastic media
Authors Jing Zeng, Alexey Stovas and Handong HuangABSTRACTAn important cause of seismic anisotropic attenuation is the interbedding of thin viscoelastic layers. However, much less attention has been devoted to layer‐induced anisotropic attenuation. Here, we derive a group of unified weighted average forms for effective attenuation from a binary isotropic, transversely isotropic‐ and orthorhombic‐layered medium in the zero‐frequency limit by using the Backus averaging/upscaling method and analyse the influence of interval parameters on effective attenuation. Besides the corresponding interval attenuation and the real part of stiffness, the contrast in the real part of the complex stiffness is also a key factor influencing effective attenuation. A simple linear approximation can be obtained to calculate effective attenuation if the contrast in the real part of stiffness is very small. In a viscoelastic medium, attenuation anisotropy and velocity anisotropy may have different orientations of symmetry planes, and the symmetry class of the former is not lower than that of the latter. We define a group of more general attenuation‐anisotropy parameters to characterize not only the anisotropic attenuation with different symmetry classes from the anisotropic velocity but also the elastic case. Numerical tests reveal the influence of interval attenuation anisotropy, interval velocity anisotropy and the contrast in the real part of stiffness on effective attenuation anisotropy. Types of effective attenuation anisotropy for interval orthorhombic attenuation and interval transversely isotropic attenuation with a vertical symmetry (vertical transversely isotropic attenuation) are controlled only by the interval attenuation anisotropy. A type of effective attenuation anisotropy for interval TI attenuation with a horizontal symmetry (horizontal transversely isotropic attenuation) is controlled by the interval attenuation anisotropy and the contrast in the real part of stiffness. The type of effective attenuation anisotropy for interval isotropic attenuation is controlled by all three factors. The magnitude of effective attenuation anisotropy is positively correlated with the contrast in the real part of the stiffness. Effective attenuation even in isotropic layers with identical isotropic attenuation is anisotropic if the contrast in the real part of stiffness is non‐zero. In addition, if the contrast in the real part of stiffness is very small, a simple linear approximation also can be performed to calculate effective attenuation‐anisotropy parameters for interval anisotropic attenuation.
-
-
-
Effects of petrophysical parameters on attenuation and dispersion of seismic waves in the simplified poroelastic theory: Comparison between the Biot's theory and viscosity‐extended theory
Authors Haixia Zhao, Bangyu Wu and Jinghuai GaoABSTRACTThe simplified macro‐equations of porous elastic media are presented based on Hickey's theory upon ignoring effects of thermomechanical coupling and fluctuations of porosity and density induced by passing waves. The macro‐equations with definite physical parameters predict two types of compressional waves (P wave) and two types of shear waves (S wave). The first types of P and S waves, similar to the fast P wave and S wave in Biot's theory, propagate with fast velocity and have relatively weak dispersion and attenuation, while the second types of waves behave as diffusive modes due to their distinct dispersion and strong attenuation. The second S wave resulting from the bulk and shear viscous loss within pore fluid is slower than the second P wave but with strong attenuation at lower frequencies. Based on the simplified porous elastic equations, the effects of petrophysical parameters (permeability, porosity, coupling density and fluid viscosity) on the velocity dispersion and attenuation of P and S waves are studied in brine‐saturated sandstone compared with the results of Biot's theory. The results show that the dispersion and attenuation of P waves in simplified theory are stronger than those of Biot's theory and appear at slightly lower frequencies because of the existence of bulk and shear viscous loss within pore fluid. The properties of the first S wave are almost consistent with the S wave in Biot's theory, while the second S wave not included in Biot's theory even dies off around its source due to its extremely strong attenuation. The permeability and porosity have an obvious impact on the velocity dispersion and attenuation of both P and S waves. Higher permeabilities make the peaks of attenuation shift towards lower frequencies. Higher porosities correspond to higher dispersion and attenuation. Moreover, the inertial coupling between fluid and solid induces weak velocity dispersion and attenuation of both P and S waves at higher frequencies, whereas the fluid viscosity dominates the dispersion and attenuation in a macroscopic porous medium. Besides, the heavy oil sand is used to investigate the influence of high viscous fluid on the dispersion and attenuation of both P and S waves. The dispersion and attenuation in heavy oil sand are stronger than those in brine‐saturated sandstone due to the considerable shear viscosity of heavy oil. Seismic properties are strongly influenced by the fluid viscosity; thus, viscosity should be included in fluid properties to explain solid–fluid combination behaviour properly.
-
-
-
Seismic dispersion and attenuation in layered porous rocks with fractures of varying orientations
Authors Denghui Xu, Tongcheng Han and Li‐Yun FuABSTRACTWave‐induced fluid flow plays an important role in affecting the seismic dispersion and attenuation of fractured porous rocks. While numerous theoretical models have been proposed for the seismic dispersion and attenuation in fractured porous rocks, most of them neglect the wave‐induced fluid flow resulting from the background anisotropy (e.g. the interlayer fluid flow between different layers) that can be normal in real reservoirs. Here, according to the theories of poroelasticity, we present an approach to study the frequency‐dependent seismic properties of more realistic and complicated rocks, i.e. horizontally and periodically layered porous rock with horizontal and randomly orienting fractures, respectively, distributed in one of the two periodical layers. The approach accounts for the dual effects of the wave‐induced fluid flow between the fractures and the background pores and between different layers (the interlayer fluid flow). Because C33 (i.e., the modulus of the normally incident P‐wave) is directly related to the P‐wave velocity widely measured in the seismic exploration, and its comprehensive dispersion and attenuation are found to be most significant, we study mainly the effects of fracture properties and the stiffness contrast between the different layers on the seismic dispersion and attenuation of C33. The results show that the increasing stiffness contrast enhances the interlayer fluid flow of the layered porous rocks with both horizontal and randomly orienting fractures and weakens the wave‐induced fluid flow between the fractures and the background pores, especially for the layered porous rock with horizontal fractures. The modelling results also demonstrate that for the considered rock construction, the increasing fracture density reduces the interlayer fluid flow while improves the dispersion and attenuation in the fracture‐relevant frequency band. Increasing fracture aspect ratio is found to reduce the dispersion and attenuation in the fracture‐relevant frequency band only, especially for the layered porous rock with horizontal fractures.
-
-
-
Seismic time–frequency analysis as a robust method to estimate the fluid saturation: A case study of carbonates reservoir, Iran
Authors Javad Jamali and Abdolrahim JavaherianABSTRACTThe propagation of seismic waves through a saturated reservoir compresses the fluid in the pore spaces. During this transition, parts of seismic energy would be attenuated because of intrinsic absorption. Rock physics models make the bridge between the seismic properties and petrophysical reality in the earth. Attenuation is one of the significant seismic attributes used to describe the fluid behaviour in the reservoirs. We examined the core samples using ultrasonic experiments at the reservoir conditions. Given the rock properties of the carbonate reservoir and experiment results, the patchy saturation mechanism was solved for substituted fluid using the theory of modulus frequency. The extracted relationship between the seismic attenuation and water saturation was used in time–frequency analysis. We performed the peak frequency method to estimate the Q factor in the Gabor domain and determined the water saturation based on the computed rock physics model. The results showed how the probable fault in the reservoir has stopped the fluid movement in the reservoir and caused touching the water‐bearing zone through drilling.
-
-
-
Dependence of borehole shear‐horizontal‐wave seismoelectric response on soil textures
Authors Weihao Wang, Zhiwen Cui, Jinxia Liu and Tribikram KunduABSTRACTThe phenomenon of acoustic waves inducing electric fields in porous media is called the seismoelectric effect. Earlier investigators proposed the usage of seismoelectric effect for well logging. Soil texture has a strong influence on the coupled wave fields during shallow surface explorations. In this article, we study the borehole pure shear‐horizontal wave and the coupling transverse‐electric field (acoustic–electrical coupling wave fields) in the partially saturated soil. Combined with related theories, we expand the formation parameters to partially saturated forms and discuss the influence of soil texture conditions on the seismoelectric wave fields. The results show that the elastic and electrical properties of porous media are sensitive to water saturation. The compositions of the acoustic and electric fields for different soil textures do not change, but the waveforms differ. We also use the secant integral method to simulate the interface‐converted electromagnetic waves. The results show that interface response strength is greatly influenced by soil texture. In addition, considering the sensitivity of the inducing electric field to fluid salinity, we also simulate the time‐domain waveforms of electric field for different pore fluid salinity levels. The results show that as the salinity increases, the electric field amplitude decreases monotonically. The above conclusions have certain significance for the application of borehole shear wave and its coupled electric fields for resource exploration, saturation assessment and groundwater pollution monitoring.
-
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)