Geophysical Prospecting - Volume 73, Issue 2, 2025
Volume 73, Issue 2, 2025
- ISSUE INFORMATION
-
- CORRIGENDUM
-
- ORIGINAL ARTICLE
-
-
-
Three‐dimensional gravity forward modelling based on rectilinear grid and Block–Toeplitz Toeplitz–Block methods
More LessAuthors Ling Wan, Shihe Li, Rui Ye, Yifan Wang, Zenghan Ma and Tingting LinABSTRACTThe main method for calculating the gravity field involves discretizing the density sources into a stack of rectangular prisms with a regular grid distribution. The analytical formulation of the gravity anomaly for a right‐angled rectangular prism is affected by depth, with the kernel function decaying as depth increases. In addition, the efficiency of the computation and the storage requirements often pose challenges. We present a fast computational method for three‐dimensional gravity forward modelling of subsurface space using rectilinear grid and apply the Block–Toeplitz Toeplitz–Block method to the rectilinear grid. The size of the upright rectangles increases with depth to offset the effect of depth. We assume that the observation points are distributed on a homogeneous grid, and the kernel sensitivity matrices exhibit a Block–Toeplitz Toeplitz–Block structure, which is symmetric. For rectilinear dissections of subsurface space in MATLAB, the logarithmic interval size is used. The rectilinear mesh can offset the effect of depth to some degree allowing gravity anomalies to decrease more quickly. For the test of a single model, the gravity anomalies decrease faster and more rapidly in the case of the rectilinear grid compared to the uniform grid. In addition to this, we performed simulations on more complex models and demonstrated that using the Block–Toeplitz Toeplitz–Block method on this basis greatly improves the computational efficiency.
-
-
-
-
Unsupervised learning inversion of seismic velocity models based on a multi‐scale strategy
More LessAuthors Senlin Yang, Bin Liu, Yuxiao Ren and Peng JiangAbstractDeep learning‐based methods have performed well in seismic waveform inversion tasks in recent years, while the need for velocity models as labels has somewhat limited their application. Unsupervised learning allows us to train the neural network without labels. When inverting seismic velocity models from observed data, labels are often unavailable for real data. To address this problem and improve network generalization, we introduce a multi‐scale strategy to enhance the performance of unsupervised learning. The first ‘multi‐scale’ is derived from the conventional full waveform inversion strategy, in which the low‐, middle‐ and high‐frequency inversion results are successively predicted during the network training. Another ‘multi‐scale’ is to introduce multi‐scale similarity as an additional data loss term to improve the inversion results. With 12,000 samples from the overthrust model, our method obtains comparable results with the supervised learning method and outperforms unsupervised methods that rely only on the mean square error as a loss function. We compare the performance of the proposed method with multi‐scale full waveform inversion on the Marmousi model, and the proposed method achieves better results at low‐ and middle‐frequencies, and, as a result, it provides good initial models for further full waveform inversion updates.
-
-
-
Approximation for reflection and transmission at a thin isotropic poroelastic bed between two isotropic elastic half‐spaces
More LessAuthors Chun Yang, Alexey Stovas, Yun Wang, Jie Qu and Xueming HeAbstractSeismic responses from a horizontal poroelastic layer provide chances to detect fluids and characterize reservoirs. The poroelastic layer can be considered a thin poroelastic bed if the layer's thickness is less than about one‐eighth of the P‐wave wavelength. Most previous theoretical studies on the reflection and transmission of waves in a model containing a thin poroelastic bed employ fluid or poroelastic medium as the overlying media. Existing approximate formulas of PP‐wave reflection coefficients are given for P‐wave normal‐incidence. Thus, this paper derived the wave reflection and transmission approximate formulas of a thin poroelastic bed between two elastic half‐spaces with P‐wave oblique incident. First, we illustrated the exact reflection and transmission matrix equations for P‐wave incidents based on poroelasticity theory and the boundary conditions. Assuming the poroelastic bed's thickness is far less than wavelengths of S‐ and P‐waves, approximate reflection and transmission formulas are expanded in Taylor series centred at value of the parameter defined as the product of angular frequency, thickness and slowness. Numerical results show that the thinner the poroelastic layer, the closer the approximate reflection and transmission coefficients are to the exact ones. The approximate formulas are valid for small and medium angles. Approximated PP‐wave reflection and transmission coefficients are closer to the exact values than those of the converted waves, which is caused by the fact that P‐wave has a lower slowness than S‐wave.
-
-
-
A workflow combining Kirchhoff single‐stack redatuming and common‐reflection surface‐based enhancement for depth imaging: A case study
More LessAuthors German Garabito, Bruno F. Gonçalves, João L. Caldeira and Hashem ShahsavaniAbstractSeismic datasets acquired in onshore basins typically present a great challenge for seismic depth imaging due to the poor quality of the prestack data caused by inherent acquisition difficulties and, in particular, the distortion of the seismic signal caused by topography and heterogeneity of the weathering zone. In this seismic data, the standard static correction does not give satisfactory results in depth imaging. Wave‐equation‐based redatuming is an alternative solution to this problem, as it correctly moves the data measured on the terrain surface to a new flat datum. But most redatuming techniques require knowledge of an accurate near‐surface velocity model. A new workflow for depth imaging of land data is proposed by combining Kirchhoff single‐stack redatuming and prestack data enhancement by common‐reflection surface stack method. A major advantage of this redatuming method is that it only requires a good approximation of the near‐surface velocity model to properly remove distortions in the seismic signal caused by the rugged topography and weathering zone. A new algorithm is introduced to apply Kirchhoff single‐stack redatuming to multi‐coverage land seismic data. Common‐reflection surface‐based prestack data enhancement is applied after redatuming to attenuate radon noise and enhance coherent event. The proposed workflow has been successfully applied to land seismic data from the Parnaíba Basin in northeastern Brazil, transforming surface‐referenced prestack data to a new flat datum. By applying common‐reflection surface‐based enhancement to the redatumed prestack data, a significant improvement in signal‐to‐noise ratio and enhancement of reflections was achieved. The depth‐migrated image confirms the great improvement in quality, where the reflectors show strong enhancement and better continuity throughout the section, compared to the migrated image obtained from the only redatumed data. This improvement in quality was important for the interpretation of the main reflectors of the geological formations.
-
-
-
Applications of deep learning‐based resolution‐enhanced seismic data in fault identification
More LessAuthors Lei Lin, Chenglong Li, Yanbin Kuang, Xing Xin and Zhi ZhongAbstractHigh‐quality seismic data play a crucial role in accurately interpreting tectonic and lithologic features such as small faults, river margins and thin beds. Over the past decades, researchers have developed numerous methods to enhance seismic resolution and signal‐to‐noise ratio. However, the benefits of quality‐improved seismic data for seismic interpretation have received limited attention. In response, we propose a generative adversarial network–based algorithm to enhance seismic quality and assess how this algorithm improves the accuracy of both machine learning–based and manual fault identification. For machine learning–based fault identification, we integrate a resolution enhancement and noise attenuation neural network (HRNet) with a fault identification neural network (FaultNet). A raw seismic image is first fed into the trained HRNet to obtain a resolution‐enhanced and noise‐suppressed image, which is then input into the trained FaultNet to produce the high‐resolution fault probability map. For manual fault identification, we enlisted three interpreters with geophysical backgrounds to annotate faults on seismic images both before and after HRNet enhancement. Comparison experiments on three field seismic samples show that our method generates more accurate, cleaner and sharper fault probability maps than directly feeding raw seismic images into FaultNet. In addition, our workflow outperforms both the milestone fault identification method and state‐of‐the‐art Transformer‐based neural networks, particularly in detecting small‐scale faults. Furthermore, the HRNet‐enhanced seismic images help interpreters identify small‐ and medium‐scale faults with reduced uncertainty. In the future, HRNet‐enhanced seismic data can be applied to a broader range of high‐precision seismic interpretation tasks, including horizon picking, channel boundary detection and attribute inversion.
-
-
-
Three‐dimensional electromagnetic inversion of transfer function data from controlled sources
More LessAuthors Paula Rulff, Thomas Kalscheuer, Mehrdad Bastani and Dominik ZbindenAbstractWe develop a three‐dimensional inversion code to image the resistivity distribution of the subsurface from frequency‐domain controlled‐source electromagnetic data. Controlled‐source electromagnetic investigations play an important role in many different geophysical prospecting applications. To evaluate controlled‐source electromagnetic data collected with complex measurement setups, advanced three‐dimensional modelling and inversion tools are required.
We adopt a preconditioned non‐linear conjugate gradient algorithm to enable three‐dimensional inversion of impedance tensor and vertical magnetic transfer function data produced by multiple sets of two independent active sources. Forward simulations are performed with a finite‐element solver. Increased sensitivities at source locations can optionally be counteracted with a weighting function in the regularization term to reduce source‐related anomalies in the resistivity model. We investigate the capabilities of the inversion code using one synthetic and one field example. The results demonstrate that we can produce reliable subsurface models, although data sets from single pairs of independent sources remain challenging.
-
-
-
Facies‐constrained simultaneous inversion for elastic parameters and fracture weaknesses using azimuthal partially incidence‐angle‐stacked seismic data
More LessAuthors Huaizhen Chen, Jian Han and Kun LiAbstractIn order to improve the identification and characterization of underground fractured reservoirs, seismic inversion for elastic properties and fracture indicators is required. To improve the accuracy of seismic inversion, model constraints are necessary. Model constraints of P‐ and S‐wave moduli can be provided by well logging data; however, model constraints of fracture weaknesses are often unavailable. To obtain model constraints of fracture weaknesses, we propose a two‐stage inversion method, which is implemented as (1) estimating azimuthal elastic impedance (AEI) and fracture facies using partially incidence‐angle‐stacked seismic data at different azimuths; and (2) using the estimated azimuthal elastic impedance to predict P‐ and S‐wave moduli, density and fracture weaknesses, which is constrained by models constructed using the estimated fracture facies. In the first stage, we use Gaussian mixture prior distribution to obtain azimuthal elastic impedance of different incidence angles and azimuths, and we also predict fracture facies combining the obtained azimuthal elastic impedance and seismic data. In the second stage, we implement the Bayesian maximum a posterior inversion for estimating unknown parameter vectors. We apply the proposed inversion method to noisy synthetic seismic data, which illustrates the inversion method is robust even in the case of a signal‐to‐noise ratio of 1. Tests on real data reveal that reliable results of P‐ and S‐wave moduli and fracture weaknesses are obtained, which verifies that the inversion method is a valuable tool for generating reliable fracture indicators from azimuthal seismic data for identifying underground fractured reservoirs.
-
-
-
A hybrid network for three‐dimensional seismic fault segmentation based on nested residual attention and self‐attention mechanism
More LessAuthors Qifeng Sun, Hui Jiang, Qizhen Du and Faming GongAbstractFault detection is a crucial step in seismotectonic interpretation and oil–gas exploration. In recent years, deep learning has gradually proven to be an effective approach for detecting faults. Due to complex geological structures and seismic noise, detection results of such approaches remain unsatisfactory. In this study, we propose a hybrid network (NRA‐SANet) that integrates a self‐attention mechanism into a nested residual attention network for a three‐dimensional seismic fault segmentation task. In NRA‐SANet, the nested residual coding structure is designed to fuse multi‐scale fault features, which can fully mine fine‐grained fault information. The two‐head self‐attention decoding structure is designed to construct long‐distance fault dependencies from different feature representation subspaces, which can enhance the understanding of the model regarding the global fault distribution. In order to suppress the interference of seismic noise, we propose a fault‐attention module and embed it into the model. It utilizes the weighted and the separate‐and‐reconstruct channel strategy to improve the model sensitivity to fault areas. Experiments demonstrate that NRA‐SANet exhibits strong noise robustness, while it can also detect more continuous and more small‐scale faults than other approaches on field seismic data. This study provides a new idea to promote the development of seismic interpretation.
-
-
-
Topography‐dependent eikonal solver in tilted transverse isotropic media with anelliptical factorization
More LessAuthors Xiaole Zhou, Serge Sambolian, Haiqiang Lan and Stéphane OpertoAbstractIncorporating anisotropy and complex topography is necessary to perform traveltime tomography in complex land environments while being a computational challenge when traveltimes are computed with finite‐difference eikonal solvers. Previous studies have taken this challenge by computing traveltimes in transverse isotropic media involving complex topography with a finite‐difference eikonal equation solver on a curvilinear grid. In this approach, the source singularity, which is a major issue in eikonal solvers, is managed with the elliptical multiplicative factorization method, where the total traveltime field is decomposed into an elliptical base traveltime map, which has a known analytical expression and an unknown perturbation field. However, the group velocity curve can deviate significantly from an ellipse in anellipitically anisotropic media. In this case, the elliptical base traveltime field differs significantly from the anelliptical counterpart, leading to potentially suboptimal traveltime solutions, even though it helps to mitigate the detrimental effects of the source singularity. To overcome this issue, we develop a more accurate topography‐dependent eikonal solver in transverse isotropic media that relies on anelliptical factorization. To achieve this, we first define the coordinate transform from the Cartesian to the curvilinear coordinate system, which provides the necessary framework to implement the topography‐dependent transverse isotropic finite‐difference eikonal solver with arbitrary source and receiver positioning. Then, we develop a semi‐analytical method for the computation of the topography‐dependent anelliptical base traveltime field. Finally, we efficiently solve the resulting quadratic elliptical equation using the fast sweeping method and a quartic anelliptical source term through fixed‐point iteration. We assess the computational efficiency, stability and accuracy of the new eikonal solver against the solver based on elliptical factorization using several transverse isotropic numerical examples. We conclude that this new solver provides a versatile and accurate forward engine for traveltime tomography in complex geological environments such as foothills and thrust belts. It can also be used in marine environments involving complex bathymetry when tomography is applied to redatumed data on the sea bottom.
-
-
-
Removing random noise and improving the resolution of seismic data using deep‐learning transformers
More LessAuthors Qifeng Sun, Yali Feng, Qizhen Du and Faming GongAbstractPost‐stack data are susceptible to noise interference and have low resolution, which impacts the accuracy and efficiency of subsequent seismic data interpretation. To address this issue, we propose a deep learning approach called Seis‐SUnet, which achieves simultaneous random noise suppression and super‐resolution reconstruction of seismic data. First, the Conv‐Swin‐Block is designed to utilize ordinary convolution and Swin transformer to capture the long‐distance dependencies in the spatial location of seismic data, enabling the network to comprehensively comprehend the overall structure of seismic data. Second, to address the problem of weakening the effective signal during network mapping, we use a hybrid training strategy of L1 loss, edge loss and multi‐scale structural similarity loss. The edge loss function directs the network training to focus more on the high‐frequency information at the edges of seismic data by amplifying the weight. Additionally, the verification of synthetic and field seismic datasets confirms that Seis‐SUnet can effectively improve the signal‐to‐noise ratio and resolution of seismic data. By comparing it with traditional methods and two deep learning reconstruction methods, experimental results demonstrate that Seis‐SUnet excels in removing random noise, preserving the continuity of rock layers and maintaining faults as well as being strong robustness.
-
-
-
Attention mechanism‐assisted recurrent neural network for well log lithology classification
More LessAuthors Yining Gao, Miao Tian, Dario Grana, Zhaohui Xu and Huaimin XuAbstractLithology classification is a fundamental aspect of reservoir classification. Due to the limited availability of core samples, computational modelling methods for lithology classification based on indirect measurements are required. The main challenge for standard clustering methods is the complex vertical dependency of sedimentological sequences as well as the spatial coupling of well logs. Machine learning methods, such as recurrent neural networks, long short‐term memory and bidirectional long short‐term memory, can account for the spatial correlation of the measured data and the predicted model. Based on these developments, we propose a novel approach using two distinct models: a self‐attention‐assisted bidirectional long short‐term memory model and a multi‐head attention‐based bidirectional long short‐term memory model. These models consider spatial continuity and adaptively adjust the weight in each step to improve the classification using the attention mechanism. The proposed method is tested on a set of real well logs with limited training data obtained from core samples. The prediction results from the proposed models and the benchmark one are compared in terms of the accuracy of lithology classification. Additionally, the weight matrices from both attention mechanisms are visualized to elucidate the correlations between depth steps and to help analyse how these mechanisms contribute to improved prediction accuracy. The study shows that the proposed multi‐head attention‐based bidirectional long short‐term memory model improves classification, especially for thin layers.
-
-
-
Elastic least‐squares reverse time migration from topography through anisotropic tensorial elastodynamics
More LessAuthors Tugrul Konuk and Jeffrey ShraggeAbstractLeast‐squares reverse time migration is an increasingly popular technique for subsurface imaging, especially in the presence of complex geological structures. However, elastic least‐squares reverse time migration algorithms face significant practical and numerical challenges when migrating multi‐component seismic data acquired from irregular topography. Many associated issues can be avoided by abandoning the Cartesian coordinate system and migrating the data to a generalized topographic coordinate system conformal to surface topology. We introduce a generalized anisotropic elastic least‐squares reverse time migration methodology that uses the numerical solutions of tensorial elastodynamics for propagating wavefields in computational domains influenced by free‐surface topography. We define a coordinate mapping assuming unstretched vertically translated meshes that transform an irregular physical domain to a regular computational domain on which calculating numerical elastodynamics solutions is straightforward. This allows us to obtain numerical solutions of forward and adjoint elastodynamics and generate subsurface images directly in topographic coordinates using a tensorial energy‐norm imaging condition. Numerical examples demonstrate that the proposed generalized elastic least‐squares reverse time migration algorithm is suitable for generating high‐quality images with reduced artefacts and better balanced reflectivity that can accurately explain observed data acquired from topography in a medium characterized by arbitrary heterogeneity and anisotropy. Finally, the computational cost of our method is comparable to that of an equivalent Cartesian elastic least‐squares reverse time migration numerical implementation.
-
-
-
Seismic imaging of deep‐seated gold deposit and host rocks through a reappraisal of legacy seismic data in the Fochville mining area, South Africa
More LessAbstractReappraisal of legacy reflection seismic data has shown to deliver value in mineral exploration, particularly in brownfield settings. In this work, we demonstrate how the reappraisal and processing of legacy reflection seismic data can be advantageous in the mineral exploration industry. We use today's standard seismic processing tools to improve the imaging of deep and complex geological structures that host mineral deposits. The recovered and processed 25.3 km long legacy seismic profile in this study was acquired in 1983 by the Gold Division of Anglo‐American as part of the Witwatersrand Gold Fields exploration program. This study aims to improve the imaging of the Ventersdorp Contact Reef gold‐bearing horizon (termed reef), a world‐class gold deposit (2 m thick) situated at depths between ∼2400 and ∼4100 m below the ground surface near South Deep Gold Mine in Fochville, South Africa. The final processing results from the pre‐stack time and phase‐shift migration approaches clearly reveal a dipping reflection associated with the gold‐bearing horizon and major steeply dipping faults that crosscut and displace the deposit. The final results are integrated with borehole information, 1D synthetic modelling and aeromagnetic data to constrain the structural interpretation. In particular, 1D synthetic simulation and borehole data constrain the depth position of the gold deposit. The magnetic data provides additional constraints on the complex faulted blocks of the host rocks such as the intrusions that may have a direct impact on ore resources and evaluation. The mining companies, such as South Deep Gold Mine, operating closer to the seismic profiles can use this new structural information to update the current geological models and improve future mine planning and designs, thus providing some insight into the prospectivity of unmined ground.
-
-
-
Longitudinal electrical resistivity tomography interpretation: Numerical modelling
More LessAuthors Balgaisha Mukanova, Marzhan Turarova and Dilyara RakishevaAbstractThis study conducts numerical simulations using the boundary integral equations method to model electric resistivity tomography for investigating an embankment dam. Computational speed is enhanced by formulating the direct problem in integral equations and employing the Fourier transform along the strike direction, thereby reducing the 2.5D problem to a set of 1D integral equations. Synthetic data obtained for potential installations are interpreted using the inversion programme ZondRes2D. The impact of varying dam slope angles and the resistivity of interfacing media on interpretation outcomes is systematically evaluated across different dam models. Main geometric anomalies associated with the three‐dimensional structure of the dam are described.
-
-
-
Deep‐learning‐based Q model building for high‐resolution imaging
More LessAuthors Xin Ju, Jincheng Xu and Jianfeng ZhangAbstractBuilding a macro Q model for deabsorption migration using surface reflection data is challenging owing to interferences of the reflections resulting from stacked thin layers. The effective Q approach gives an alternative way to overcome this difficulty. However, manual processing is involved for effective Q estimation. This restricts the use of denser grids in building an inhomogeneous Q model. We therefore incorporate deep learning into the effective Q approach, thus yielding a deep learning‐based Q model building scheme. The resulting scheme improves the manual effective Q estimation by simultaneously accounting for the imaging resolution and induced noises using two networks. Moreover, most manual processing is reduced in spite of denser grids in building a 3D Q model. One of the networks used is a 1D convolutional neural network that determines the optimal upper cut‐off frequency for a selected Q with an input of multi‐channel amplitude spectra, and another is a residual neural network that determines the optimal Q for a series of Q values with an input of multi‐channel imaging sections inside the selected small window filtered under the corresponding upper cut‐off frequencies. As a result, a Q model that improves the imaging resolution in the absence of amplification of noises is gained. Transfer learning is used, thus reducing the training cost when applied to different geological targets. We test our scheme using 3D field data. Higher resolution images without induced noises are obtained by a deabsorption migration using the Q model built and compared to those obtained by the migration without absorption compensation.
-
-
-
A dual‐porosity model incorporating uniaxial stress effect and its application in wave velocity estimation with well‐logging data
More LessAuthors Jiayun Li, Zhaoyun Zong and Fubin ChenAbstractThe elastic properties and wave propagation of porous rocks are sensitive to the stress variation. The existing theories mainly focus on the impacts of effective stress, confining and pore pressures. The physics of uniaxial stress effect on rock elasticity and wave propagation is seldom well studied, although the uniaxial stress case is frequently encountered in several scenarios, such as the laboratory loading and the subsurface tectonic deformation. Therefore, we propose a new dual‐porosity model to describe the effect of uniaxial effective stress on elastic properties of dry porous rocks, based on the Palmer equation and Shapiro dual‐porosity model. The Gurevich squirt‐flow model is then incorporated to model the dispersion and attenuation of wave velocities of fluid‐saturated porous rocks. Modelling results show that the increase of uniaxial effective stress inflates the P‐ and S‐wave velocities along the stress direction until the velocities asymptotically reach their maximum values within the elastic limit. However, the relevant wave dispersion and attenuation gradually decline with the elevating stress possibly due to the gradual closure of cracks. The effect of viscosity, fluid modulus and crack aspect ratio on wave dispersion is investigated in detail as well. By comparing our model to the published laboratory ultrasonic measurements, we confirm the validity of our model. Furthermore, our dual‐porosity model is used to establish a rock‐physics approach to estimate the wave velocities with the well‐logging data. The well‐logging examples show the reasonable agreement between the predicted results and real data, illustrating the feasibility of our approach.
-
Volumes & issues
-
Volume 73 (2024 - 2025)
-
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)
Most Read This Month