Exploration Geophysics - Volume 49, Issue 3, 2018
Volume 49, Issue 3, 2018
- Reasearch Articles
-
-
-
Noise reduction of grounded electrical source airborne transient electromagnetic data using an exponential fitting-adaptive Kalman filter
More LessAuthors Yanju Ji, Qiong Wu, Yuan Wang, Jun Lin, Dongsheng Li, Shangyu Du, Shengbao Yu and Shanshan Guan[GREATEM field data usually includes a mixed variety of noises, which makes the exponential decaying signal too difficult to identify. This paper presents an exponential fitting-adaptive Kalman filter (EF-AKF) to remove mixed electromagnetic noises, while preserving the signal characteristics. As a new method, the EF-AKF can be used for denoising exponential decaying signals.
,The grounded electrical source airborne transient electromagnetic (GREATEM) system, which uses a grounded electrical transmitter and an aircraft for the receiver, offers deep exploration capability and detection efficiency. However, GREATEM field data usually includes mixed varied noises (white noise, sferics noise and human noise), which make identifying the exponential decaying signal too difficult. Traditional filtering methods mainly focus on suppressing specific noise types, which may cause the distortion of GREATEM signal, especially when the signal is affected by high residual sferics noise. This paper presents an exponential fitting-adaptive Kalman filter (EF-AKF) to remove mixed electromagnetic noises, while preserving the signal characteristics. The EF-AKF consists of an exponential fitting procedure and an adaptive scalar Kalman filter (SKF). The adaptive SKF uses the exponential fitting results in the weighting coefficients calculation. The EF-AKF is verified on an analytical three-layer model. It is compared with the SKF and wavelet threshold-exponential adaptive window width-fitting denoising algorithm (WEF) in synthetic data. The results showed that the EF-AKF outperformed the other methods in the noise reduction of GREATEM data. The EF-AKF is also tested on a synthetic quasi-2D earth model and applied to GREATEM field data in Huaide, Jilin province, China. Application of the EF-AKF allowed considerable improvement of the quality of the GREATEM field data.
]
-
-
-
-
Three-dimensional inversion of CSAMT data in the presence of topography
More LessAuthors Changhong Lin, Handong Tan, Wangwang Wang, Tuo Tong, Miao Peng, Mao Wang and Weihua Zeng[In this paper, we present a scheme to incorporate 3D controlled-source audio frequency magnetotelluric (CSAMT) topographic distortions into the 3D inversion instead of correcting them. This approach has been verified by comparison with 2D FEM CSAMT solutions and synthetic inversion examples. The field example also illustrates the effectiveness of our approach.
,3D controlled-source audio frequency magnetotelluric (CSAMT) responses can be distorted strongly by topography and should be accounted for in data inversion and interpretation. In this paper we present a scheme to incorporate topographic distortions into the inversion instead of correcting them. This approach has been verified by comparing the modelling results with 2D FEM CSAMT solutions and synthetic inversion examples. Compared with the responses generated from a half-space model with flat surface, it is found that not only the topography in the survey area but also that at the source position may strongly distort the CSAMT responses. The field example indicates that results with topography are much better than those without considering topography to map the distribution of coal seam underground, which also illustrates the effectiveness of our approach.
]
-
-
-
Three-dimensional tensor controlled-source audio-frequency magnetotelluric inversion using LBFGS
More LessAuthors Kunpeng Wang, Handong Tan, Changhong Lin, Jianlong Yuan, Cong Wang and Jing Tang[We have shown in this paper that directly using a magnetotelluric method to invert the data of the tensor controlled-source audio-frequency magnetotelluric (CSAMT) will obtain an incorrect result, the inversion result of tensor CSAMT is more reliable than that of the traditional CSAMT, and the limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method is more efficient than the nonlinear conjugate gradient (NLCG) method for tensor CSAMT.
,The controlled-source audio-frequency magnetotelluric (CSAMT) method has become an important method in geophysical electromagnetic exploration. However, traditional CSAMT only gathers a single set of orthogonal electric and magnetic data, which cannot describe the whole subsurface geological structure. Due to increasingly complex geological targets, the drawbacks of traditional CSAMT have gradually become more significant, promoting the need for tensor CSAMT. Tensor CSAMT can gather richer information, but the 3D forward and inversion models of this method have developed slowly since it was first proposed. The common method for inverting the data of the tensor CSAMT is still magnetotelluric (MT). This paper adopts a staggered-grid finite difference method to realise the 3D forward modelling of the tensor CSAMT. On this basis, we adopt a limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) method to implement a 3D inversion with full impedance data. Through inverting synthetic and real data, we prove that: (1) directly using an MT method to invert the data of the tensor CSAMT will obtain an incorrect result, (2) the inversion result of tensor CSAMT is more reliable than that of the traditional CSAMT, and (3) LBFGS is more efficient than the nonlinear conjugate gradient (NLCG) for tensor CSAMT. Our research shows that 3D tensor CSAMT inversion with LBFGS is very useful and practical for electromagnetic exploration.
]
-
-
-
Mapping electrical structures in the southern Great Khingan Range, north-east China, through two-dimensional magnetotelluric sounding
More LessAuthors Weijun Zhao, Shuwang Chen and Qiuhong Ding[Magnetotelluric measurements were carried out in the southern Great Khingan Range, Inner Mongolia, China, from 2012 to 2014. After 2D inversion, the Paleozoic Hunnitu and Gadasu depressions were discovered and inferred to possibly contain Linxi strata. The Mesozoic Xiretu depression was also discovered, which possibly contains Cretaceous hydrocarbon-bearing strata.
,Magnetotelluric (MT) measurements were performed in the southern Great Khingan Range from 2012 to 2014 to reveal the tectonic framework of the Jarud rift basin, and to identify the strata distribution of the Permian Linxi Formation. The dark mudstones in the Linxi Formation have shown the potential of shale gas. The study area is geographically situated in the Jarud Banner and Ar Horqin Banner, Inner Mongolia, China. It is covered predominantly with Cretaceous–Jurassic igneous rocks except for the small south-eastern part. It is tectonically located in the southern Great Khingan Range on the western periphery of the Songliao Basin, north of the Xar Moron Fault.
MT data were acquired on six north–west and five north–east profiles in the study area. After geoelectric dimensionality analysis, strike analysis and 2D inversion, the resulting resistivity models revealed numerous large faults, part of which formed the borders of the Jarud Basin. In contrast to the original extent of the Jurassic Jarud Basin, the new borders should extend further to the north. According to well logging data and geological information, two Paleozoic depressions inside the Jarud Basin, Hunnitu depression and Gadasu depression were discovered and inferred to possibly contain Linxi strata. Meanwhile, the Mesozoic Xiretu depression outside of the Jarud Basin was discovered and inferred to possibly contain Cretaceous hydrocarbon-bearing strata. In particular, the Gadasu depression, with an area of ~500 km2, contains reasonably thick conductive sediments exceeding 4 km in depth, which were inferred to be dark mudstones possibly pertaining to shale gas. Overall, the MT work will be helpful for the assessment of potential shale gas.
]
-
-
-
Estimating high hydraulic conductivity locations through a 3D simulation of water flow in soil and a resistivity survey
More Less[In this study, we propose a method to estimate high hydraulic conductivity locations that uses 3D simulation of soil water flow and in-line resistivity survey data acquired during a groundwater recharge experiment, and we apply this method to numerical and field experiments. The high hydraulic conductivity locations are estimated from a combination of field-observed and simulated apparent resistivities using the following simple steps. (1) Assuming that high hydraulic conductivity zones exist in the first layer, simulations of saturated-unsaturated seepage are conducted for several possible water-flow models that have high hydraulic conductivity zones in different locations. (2) The simulated volumetric water contents are converted into bulk resistivities, which are used to produce apparent resistivity data through simulation of a resistivity survey. (3) The differences between the simulated apparent resistivities and the field-observed data are examined, and the best-fit hydraulic conductivity model is identified by minimising the above differences. In the numerical experiment, 3D inversion of the simulated resistivity survey provides an image of the preferential flow, although the infiltration locations are unclear. Comparing the field model with the possible models, the high hydraulic conductivity location in the field model corresponds to the high hydraulic conductivity location in the possible model with the minimum errors. In the field, an in-line resistivity survey was conducted during a groundwater recharge experiment on a pyroclastic plateau. The 3D inversion of the in-line resistivity survey data provides an image of the preferential flow. Comparing the field apparent resistivity data with the simulated apparent resistivity data, the high hydraulic conductivity location of the possible model that provides the minimum error corresponds to the recharge water range, whereas the hydraulic conductivity location of the possible model that gives the maximum errors corresponds to ranges with no recharge water. These results indicate that it is possible to estimate high hydraulic conductivity locations using 3D simulations of the soil water flow and a resistivity survey.
,We propose a simple method for estimating high hydraulic conductivity locations. The proposed method uses the 3D simulations of soil water flow and resistivity survey during a groundwater recharge experiment. Results of numerical and field experiments indicate that the proposed method estimates the high hydraulic conductivity locations more precisely compared with 3D inversion of in-line data.
]
-
-
-
A novel approach to comparing AEM inversion results with borehole conductivity logs
More LessAuthors Niels B. Christensen and Kenneth C. Lawrie[Two consistent methods of comparing borehole induction logs with models from inversion of AEM data have been developed: one in model space and one in data space. The two methods are applied to the Broken Hill Managed Aquifer Recharge project conducted by Geoscience Australia on AEM data from the SkyTEM system.
,Borehole conductivity logs, besides being useful for identifying, interpreting and correlating geological formations, also find widespread use as auxiliary information in the inversion of airborne electromagnetic (AEM) data. One of the quality checks often applied to AEM inversion results is a comparison between the conductivity structures revealed by borehole conductivity logs in the survey area and the AEM inversion model closest to the borehole, often called an ‘FID point comparison’.
Another use of borehole conductivity logs is found in modern AEM inversion procedures, where the borehole conductivity information is included as prior information in a laterally constrained inversion. In most former and present practices, AEM layer conductivities are compared with the measured conductivity in the borehole. However, the borehole conductivity is essentially an apparent conductivity – it is a measured data value – while the AEM layer conductivities are model parameters resulting from inverting AEM data. To avoid comparing data and model parameters we suggest a conceptually clear approach based on an inversion of the borehole conductivity data to obtain a borehole conductivity model, which in turn can be compared with the AEM model. Furthermore, the AEM forward response of the borehole model can, in a consistent way, be compared with the AEM data. In both approaches, we keep track of uncertainty and define quantitative, uncertainty-normalised measures of the difference between borehole and AEM values, and we find simple functional relationships between the two. The methodology is demonstrated on the AEM data and conductivity logs of the Broken Hill Managed Aquifer Recharge (BHMAR) project.
]
-
-
-
Regularisation parameter adaptive selection and its application in the prestack AVO inversion
More LessAuthors Guangtan Huang, Jingye Li, Cong Luo and Xiaohong Chen[The numerical solution of the inverse problem is usually obtained by solving a set of linear algebraic equations, while the system of equations may suffer from ill-posedness due to insufficient data. Regularisation is a technique for making the estimation problems well posed by adding indirect constraints on the estimated model, but the regularisation parameter selection is difficult. In geophysics, without explicit calculation methods and quantitative evaluation criteria, it is usually based on the experience of the inversion engineers to try to achieve the best inversion results by continuously modifying the regularisation parameter. For prestack amplitude variation with offset (AVO) inversion for real seismic data, fixed regularisation parameters cannot satisfy the optimisation conditions in the seismic data with different signal-to-noise (SNR) in one area. Besides, fixed regularisation parameter may cause that the model constraint misfit is too large or too small compared to data misfit, which may guide the inversion to generate undesirable results. Therefore, adaptive selection of regularisation parameter according to the seismic data can help guarantee a good inversion result. Based on the traditional L-curve criterion, we derive the theoretical formula of the adaptive computation of the regularisation parameter, which can be applied to any norm constraint. We proposed the application of this selection scheme in prestack AVO inversion. Model tests show that the improved L-curve method has better stability than its main competitor, the generalised cross-validation (GCV) method. Prestack AVO inversion on logging data and real seismic data demonstrate that the proposed method can improve the accuracy of the inversion, and it is more immune to strong noise.
,In this paper, based on L-curve criterion, we propose an improved method for the adaptive acquisition of regularisation parameters for arbitrary norm condition. A detailed derivation of the proposed method is described. Numerical experiments confirm that the proposed method is more accurate and robust than its main competitor, generalised cross-validation.
]
-
-
-
A robust surface-consistent residual phase correction method based on migrated gathers
More LessAuthors Jincheng Xu, Hao Zhang and Jianfeng Zhang[Conventional residual static corrections determine the residual statics (time shifts) using common mid-point (CMP) gathers to estimate and correct anomalies induced by the near-surface. These time shifts disregard the phase errors existing in the data, which will reduce the resolution of the stacked image. Further errors result when reflection events in CMP gathers do not exhibit hyperbolic moveout. In order to solve these problems, we propose a robust surface-consistent residual phase correction method that simultaneously resolves both time shifts and constant phase rotations based on migrated gathers. The surface-consistent residual statics and phase are obtained from the migrated gathers expressed in terms of shot and receiver locations. We modified the standard technique of estimating the time shift corrections to include a surface-consistent constant phase rotation term. The proposed dual parameter algorithm (time shift and constant phase rotation) proved on a synthetic example that it was superior to conventional residual statics for improving the coherence of trace gathers. The computational effort can be reduced by generating migrated gathers and estimating the dual parameters in a spatially varying time window. We applied the proposed method to both synthetic and real data, and improved results were obtained with both.
,We present a novel approach for surface-consistent residual phase corrections based on migrated gathers to improve the SNR and resolution of migrated images. This method includes the following new aspects: the residual statics and phase are estimated by dual parameter cross-correlation, and determined after migration.
]
-
-
-
3D seismic attributes for structural mapping and enhancement of deep gold mining: a case study from the West Wits Line goldfields, South Africa
More LessAuthors Nomqhele Z. Nkosi, Musa S. D. Manzi, Oleg Brovko and Raymond J. Durrheim[In this study, we use sophisticated seismic attributes to interpret the 3D reflection seismic data to: (1) assess and mitigate the risks posed by deep mining activities and (2) improve the resource evaluation of the gold-bearing quartz pebble conglomerate horizons (reefs) in the world’s deepest gold mines (South Africa).
,Volumetric and horizon-based seismic attributes were carried out on the 3D reflection seismic data acquired on the deep crystalline rock environments of the Witwatersrand Basin (South Africa). This study aimed at improving the delineation of minor faults that cross-cut the gold-bearing quartz pebble conglomerate horizons (termed reefs) in the complex geology of the West Wits Line (WWL) goldfield using different seismic attributes. The 3D seismic data show delineation of faults that cross-cut and displace the target horizons by tens to hundreds of metres, leading to the delineation of the ore resources in faulted areas of the mines. In particular, the ant-tracking technique has provided better stratigraphic and structural imaging in complex faulted areas, relative to conventional interpretation methods. Other improvements include more accurate mapping of the depths, dip and strike of the key seismic horizon (Roodepoort shale) using the edge detection attributes, yielding a better understanding of the interrelationship between fault activity, reef distribution and the relative chronology of tectonic events. The integration of the 3D reflection seismic data, seismic attributes and information from borehole logs and underground mapping has provided better imaging and modelling of important fault systems that might have a direct effect on mining. This information could be used for future mining planning and designs to: (1) assess and mitigate the risks posed by mining activities, and (2) improve the resource evaluation of the gold-bearing reefs in the WWL goldfield.
]
-
-
-
Azimuthal seismic responses from shale formation based on anisotropic rock physics and reflectivity method: a case study from south-west China
More Less[We have constructed an anisotropic rock physics model at the seismic scale for shales, and applied it to the Longmaxi Shale in the Sichuan Basin, south-west China. The simplified reflectivity method is proposed to calculate seismic responses for amplitude variation with offset and azimuth (AVAz) analysis.
,Due to intrinsic anisotropy related to preferred alignment of clay particles and the existence of vertical or high angle fractures, shales usually present orthorhombic anisotropy. The objective of this study was to build anisotropic rock physics models for shales at the seismic scale. Based on the well-log and Formation Micro Imager (FMI) log data from a shale formation in the Sichuan Basin in south-west China, we derive an orthorhombic model at the seismic scale by using Schoenberg and Helbig’s method and generalised Backus averaging method in the rock physics workflow. In order to understand the relationship between physical properties of the rock physics model and seismic wave propagation, we apply the simplified reflectivity method to calculate seismic responses for amplitude variation with azimuth (AVAz) analysis. The method is based on the scheme of anisotropic reflectivity method which is commonly used to simulate the full-wave field in stratified anisotropic media, in analogy with the formula of horizontal slowness components in Schoenberg and Protázio’s method. The AVAz analysis is conducted on the seismograms of PP-wave, radial and transverse components of PS-wave. The results show that overburden effects caused by wave propagation in anisotropic media can’t be ignored. Azimuthal variations in amplitudes of both PP-wave and radial component of PS-wave can be used to indicate strikes of fractures, while PS-wave appears to be more sensitive.
]
-
-
-
Study of Scholte wave dispersion curves and modal energy distribution using a wavefield numerical simulation method
More LessAuthors Hua Wu, Guang-Zhou Shao and Qing-Chun Li[In this paper, we examine the characteristics of Scholte wave dispersion curves and their modal energy distribution by combining analytical solutions of the dispersion equation with numerical modelling. The modelling results showed that our approach can facilitate Scholte wave exploration in water-covered areas.
,Surface waves travelling along the subsurface interface between water and sediment are called Scholte (or generalised Rayleigh) waves. In this paper, we examine the characteristics of Scholte wave dispersion curves and their modal energy distribution through synthetic simulation tests. Analytical solutions of the dispersion equation for a homogeneous isotropic layered medium describe the propagation characteristics of the fundamental and higher modes, but do not provide information about the relative energy distributions of those wave modes. In real situations, the dispersion curves with different modes have different energies corresponding to different frequency bands. It is instructive to study and analyse the behaviours of Scholte wave dispersion curves by using joint dispersion equation solutions and numerical simulations results. Dispersion equation solutions are used to validate the numerical results, which in turn are used to assess energy distribution for both the fundamental and the higher modes. We observe that the thickness of the overlying water layer and the elastic properties of the ocean bottom play important roles in the dispersion characteristics of Scholte waves. The dispersion of the soft ocean-bottom model (whose S-wave velocity of the first solid layer is lower than the P-wave velocity of water) is dominated by the fundamental mode and appears to be less dependent on the thickness of the overlying water layer. In contrast, the dispersion of the hard ocean-bottom model contains significant energy in the higher modes. Higher modes will shift to lower frequencies and will be spatially closer (in the frequency direction) when the water layer thickness increases. The dispersion characteristics of the real data are in a good agreement with the conclusion of a soft ocean-bottom model. Our modelling approach of combining analytical solutions of the dispersion equation with numerical modelling can facilitate Scholte wave exploration in water-covered areas.
]
-
-
-
Acoustic wave propagation simulation by reduced order modelling
More Less[A reduced order modelling method is introduced for wave propagation modelling by using the mode shapes of the model. The numerical accuracy, computational performance and boundary conditions of the proposed method were investigated. According to the results, the proposed method substantially improves the computational efficiency of the reverse time migration.
,Wave propagation simulation, as an essential part of many algorithms in seismic exploration, is associated with high computational cost. Reduced order modelling (ROM) is a known technique in many different applications that can reduce the computational cost of simulation by employing an approximation of the model parameters. ROM can be carried out using different algorithms. The method proposed in this work is based on using the most important mode shapes of the model, which can be computed by an efficient numerical method. The numerical accuracy and computational performance of the proposed method were investigated over a simple synthetic velocity model and a portion of the SEG/EAGE velocity model. Different boundary conditions were discussed, and among them the random boundary condition had higher performance for applications like reverse time migration (RTM). The capability of the proposed method for RTM was evaluated and confirmed by the synthetic velocity model of SEG/EAGE. The results showed that the proposed ROM method, compared with the conventional finite element method (FEM), can decrease the computational cost of wave propagation modelling for applications with many simulations like the reverse time migration. Depending on the number of simulations, the proposed method can increase the computational efficiency by several orders of magnitudes.
]
-
-
-
Regional and residual gravity anomaly separation using the singular spectrum analysis-based low pass filtering: a case study from Nagpur, Maharashtra, India
More LessAuthors K. Satish Kumar, Rekapalli Rajesh and Rama Krishna Tiwari[We present here a singular spectrum analysis (SSA)-based low pass filtering algorithm for regional and residual gravity anomaly separation. The data adaptive decomposition in SSA frequency filtering algorithm facilitates reduction of the artefacts in the filtering of the non-linear and non-stationary gravity data. Initially, the method was tested on synthetic gravity data representing combined response of regional and residual components derived from the geological structures like infinite horizontal layers, intrusive dykes, deep-seated faults and volcanic intrusive bodies and compared with fast Fourier transform (FFT) wavenumber and wavelet-based filtering techniques. The results show that SSA-based filtered output exhibits a better match with pure synthetic data than the output generated from the FFT and wavelet filtering methods. The underlying method was then applied to two parallel gravity profiles of real data from the Umred coalfield, Nagpur, Maharashtra, India. Further, the SSA-based filtered regional anomaly was modelled using Geosoft GM-SYS software to construct the model of crustal structure of the study region. In essence, the modelling results suggest the following: (i) a basin-type graben structure with variable trap and sedimentary thicknesses and (ii) deep seated faults on either sides of the basin. These results correlate fairly well with known regional geology attesting the authenticity of the regional models generated from the two gravity profiles, which also agree well with each other. We, therefore, conclude that the SSA-based filtering technique is robust for regional gravity anomaly separation and could be effectively exploited for filtering other geophysical data.
,We present a singular spectrum analysis (SSA)-based algorithm for regional and residual gravity anomaly separation. Initially, we tested the method on synthetic data and applied gravity profiles from the Umred coalfield, Nagpur, Maharashtra, India. Final 2D crustal models obtained from the filtered regional anomaly agree well with regional geology and borehole information.
]
-
-
-
Enhancement of fault interpretation using multi-attribute analysis and artificial neural network (ANN) approach: a case study from Taranaki Basin, New Zealand
More LessAuthors Priyadarshi Chinmoy Kumar and Animesh Mandal[This paper describes an improved and efficient workflow that refines the art of structural interpretation from 3D seismic data. The research aims to streamline the interpretation workflow using better data conditioning scheme and defining a meta-attribute (fault probability cube) using an artificial neural network. A robust pathway is provided for delineating structural details of complex geological terrain.
,Enhanced seismic data conditioning and multi-attribute analysis through non-linear neural processing workflows has been applied to 3D seismic data over 215.10 km2 area of the Opunake prospect located in the south-eastern offshore Taranaki Basin. The present work aims to delineate faults and the related detail of structural features from the study area. Post-stack seismic data conditioning techniques such as dip-steering and structural filtering are applied to enhance the lateral continuity of seismic events and eliminate random noises from the data with the objective of improving the visibility of faults in the data volume. The conditioned data is then used to extract several attributes, such as similarity, dip variance, curvature, energy and frequency, that act as potential contributors for enhancing the fault visibility. A fully connected multilayer perceptron (MLP) network is developed to choose the proper combination of attributes for fault detection. These seismic attributes (known as the test datasets) are then trained at identified fault and non-fault locations using this network. The network comprises of 11, 7, and 2 nodes in the input layer, hidden layer and output layer, respectively. The neural training resulted in an overall minimum root mean square (RMS) misfit and misclassification (%) ranging from 0.54 to 0.67 and 18.67 to 10.42, respectively, between the trained and the test datasets. The neural training generates a fault probability attribute that produces an improved fault visibility capturing the minute details of the seismic volume as compared with the results of individual seismic attribute. Thus, the present work demonstrates an enhanced and robust workflow of fault prediction and visualisation for detail structural interpretation from 3D seismic data volume.
]
-
Volumes & issues
-
Volume 56 (2025)
-
Volume 55 (2024)
-
Volume 54 (2023)
-
Volume 53 (2022)
-
Volume 52 (2021)
-
Volume 51 (2020)
-
Volume 50 (2019)
-
Volume 49 (2018)
-
Volume 48 (2017)
-
Volume 47 (2016)
-
Volume 46 (2015)
-
Volume 45 (2014)
-
Volume 44 (2013)
-
Volume 43 (2012)
-
Volume 42 (2011)
-
Volume 41 (2010)
-
Volume 40 (2009)
-
Volume 39 (2008)
-
Volume 38 (2007)
-
Volume 37 (2006)
-
Volume 36 (2005)
-
Volume 35 (2004)
-
Volume 34 (2003)
-
Volume 33 (2002)
-
Volume 32 (2001)
-
Volume 31 (2000)
-
Volume 30 (1999)
-
Volume 29 (1998)
-
Volume 28 (1997)
-
Volume 27 (1996)
-
Volume 26 (1995)
-
Volume 25 (1994)
-
Volume 24 (1993)
-
Volume 23 (1992)
-
Volume 22 (1991)
-
Volume 21 (1990)
-
Volume 20 (1989)
-
Volume 19 (1988)
-
Volume 18 (1987)
-
Volume 17 (1986)
-
Volume 16 (1985)
-
Volume 15 (1984)
-
Volume 14 (1983)
-
Volume 13 (1982)
-
Volume 12 (1981)
-
Volume 11 (1980)
-
Volume 10 (1979)
-
Volume 9 (1978)
-
Volume 8 (1977)
-
Volume 7 (1976)
-
Volume 6 (1975)
-
Volume 5 (1974)
-
Volume 4 (1973)
-
Volume 3 (1972)
-
Volume 2 (1971)
-
Volume 1 (1970)
Most Read This Month