基于PSO的地震盲反褶积方法

本文基于地层反射系数非高斯的统计特性,在反褶积输出单位方差约束下,将反褶积输出的负熵表示为非多项式函数,作为盲反褶积的目标函数,然后采用粒子群算法优化目标函数寻找最佳反褶积算子,实现地震信号的盲反褶积.数值模拟和实际资料处理结果表明,与传统反褶积方法相比,本文方法同时适应于最小相位子波及混合相位子波的反褶积,能够更好地从地震数据中估计反射系数,有效拓宽地震资料的频谱,得到高分辨率的地震资料.     

基于带状混合矩阵ICA实现地震盲反褶积

基于对地震反褶积本质上是一个盲过程的认识,引入高阶统计学盲源分离技术——独立分量分析（ICA）实现地震盲反褶积.在无噪声假设条件下,利用地震记录时间延迟矩阵和地震子波带状褶积矩阵,将地震褶积模型转化为一般线性混合ICA模型,采用FastICA算法,将带状性质作为先验信息,实现所谓带状ICA算法（B\|ICA）,得到个数与子波算子长度相等的多个估计反射系数序列和估计子波序列,最后利用褶积模型提供的附加信息从中优选出最佳的反射系数序列及相应的地震子波.模型数据和实际二维地震道数值算例表明：对于统计性反褶积,在不对反射系数作高斯白噪假设,不对子波作最小相位假设的所谓“全盲”条件下,基于ICA方法（反射系数非高斯分布,地震子波非最小相位）可以较好解决地震盲反褶积问题,是基于二阶统计特性的地震信号统计性反褶积方法的提升,具有可行性和应用前景.     

Nonstationary deconvolution using maximum kurtosis optimization

Deconvolution is an essential step for high-resolution imaging in seismic data processing. The frequency and phase of the seismic wavelet change through time during wave propagation as a consequence of seismic absorption. Therefore, wavelet estimation is the most vital step of deconvolution, which plays the main role in seismic processing and inversion. Gabor deconvolution is an effective method to eliminate attenuation effects. Since Gabor transform does not prepare the information about the phase, minimum-phase assumption is usually supposed to estimate the phase of the wavelet. This manner does not return the optimum response where the source wavelet would be dominantly a mixed phase. We used the kurtosis maximization algorithm to estimate the phase of the wavelet. First, we removed the attenuation effect in the Gabor domain and computed the amplitude spectrum of the source wavelet; then, we rotated the seismic trace with a constant phase to reach the maximum kurtosis. This procedure was repeated in moving windows to obtain the time-varying phase changes. After that, the propagating wavelet was generated to solve the inversion problem of the convolutional model. We showed that the assumption of minimum phase does not reflect a suitable response in the case of mixed-phase wavelets. Application of this algorithm on synthetic and real data shows that subtle reflectivity information could be recovered and vertical seismic resolution is significantly improved.     

基于广义高斯分布的地震盲反褶积方法研究

Bussgang算法是针对褶积盲源分离问题提出的,本文将其用于地震盲反褶积处理.由于广义高斯概率密度函数具有逼近任意概率密度函数的能力,从反射系数序列的统计特征出发,引入广义高斯分布来体现反射系数序列超高斯分布特征.依据反射系数序列的统计特征和Bussgang算法原理,建立以Kullback-Leibler距离为非高斯性度量的目标函数,并导出算法中涉及到的无记忆非线性函数,最终实现了地震盲反褶积.模型试算和实际资料处理结果表明,该方法能较好地适应非最小相位系统,能够同时实现地震子波和反射系数估计,有效地提高地震资料分辨率.     

CORRECTING FOR COLOURED PRIMARY REFLECTIVITY IN DECONVOLUTION1

Statistical deconvolution, as it is usually applied on a routine basis, designs an operator from the trace autocorrelation to compress the wavelet which is convolved with the reflectivity sequence. Under the assumption of a white reflectivity sequence (and a minimum-delay wavelet) this simple approach is valid. However, if the reflectivity is distinctly non-white, then the deconvolution will confuse the contributions to the trace spectral shape of the wavelet and reflectivity. Given logs from a nearby well, a simple two-parameter model may be used to describe the power spectral shape of the reflection coefficients derived from the broadband synthetic. This modelling is attractive in that structure in the smoothed spectrum which is consistent with random effects is not built into the model. The two parameters are used to compute simple inverse- and forward-correcting filters, which can be applied before and after the design and implementation of the standard predictive deconvolution operators. For whitening deconvolution, application of the inverse filter prior to deconvolution is unnecessary, provided the minimum-delay version of the forward filter is used. Application of the technique to seismic data shows the correction procedure to be fast and cheap and case histories display subtle, but important, differences between the conventionally deconvolved sections and those produced by incorporating the correction procedure into the processing sequence. It is concluded that, even with a moderate amount of non-whiteness, the corrected section can show appreciably better resolution than the conventionally processed section.     

非稳态地震稀疏约束反褶积研究（英文）

传统Robinson褶积模型主要受缚于三种不合理的假设,即白噪反射系数、最小相位地震子波与稳态假设,而现代反射系数反演方法（如稀疏约束反褶积等）均在前两个假设上寻求突破的同时却忽视了一个重要事实：实际地震信号具有典型的非稳态特征,这直接冲击着反射系数反演中地震子波不随时间变化的这一基础性假设。本文首先通过实际反射系数测试证实,非稳态效应造成重要信息无法得到有效展现,且对深层影响尤为严重。为校正非稳态影响,本文从描述非稳态方面具有普适性的非稳态褶积模型出发,借助对数域的衰减曲线指导检测非稳态影响并以此实现对非稳态均衡与校正。与常规不同,本文利用对数域Gabor反褶积仅移除非稳态影响,而将分离震源子波和反射系数的任务交给具有更符合实际条件的稀疏约束反褶积处理,因此结合两种反褶积技术即可有效解决非稳态特征影响,又能避免反射系数和地震子波理想化假设的不利影响。海上地震资料的应用实际表明,校正非稳态影响有助于恢复更丰富的反射系数信息,使得与地质沉积和构造相关的细节特征得到更加清晰的展现。     

Sparsity‐enhanced wavelet deconvolution

We propose a three‐step bandwidth enhancing wavelet deconvolution process, combining linear inverse filtering and non‐linear reflectivity construction based on a sparseness assumption. The first step is conventional Wiener deconvolution. The second step consists of further spectral whitening outside the spectral bandwidth of the residual wavelet after Wiener deconvolution, i.e., the wavelet resulting from application of the Wiener deconvolution filter to the original wavelet, which usually is not a perfect spike due to band limitations of the original wavelet. We specifically propose a zero‐phase filtered sparse‐spike deconvolution as the second step to recover the reflectivity dominantly outside of the bandwidth of the residual wavelet after Wiener deconvolution. The filter applied to the sparse‐spike deconvolution result is proportional to the deviation of the amplitude spectrum of the residual wavelet from unity, i.e., it is of higher amplitude; the closer the amplitude spectrum of the residual wavelet is to zero, but of very low amplitude, the closer it is to unity. The third step consists of summation of the data from the two first steps, basically adding gradually the contribution from the sparse‐spike deconvolution result at those frequencies at which the residual wavelet after Wiener deconvolution has small amplitudes. We propose to call this technique “sparsity‐enhanced wavelet deconvolution”. We demonstrate the technique on real data with the deconvolution of the (normal‐incidence) source side sea‐surface ghost of marine towed streamer data. We also present the extension of the proposed technique to time‐varying wavelet deconvolution.     

Klauder wavelet removal before vibroseis deconvolution

The spiking deconvolution of a field seismic trace requires that the seismic wavelet on the trace be minimum phase. On a dynamite trace, the component wavelets due to the effects of recording instruments, coupling, attenuation, ghosts, reverberations and other types of multiple reflection are minimum phase. The seismic wavelet is the convolution of the component wavelets. As a result, the seismic wavelet on a dynamite trace is minimum phase and thus can be removed by spiking deconvolution. However, on a correlated vibroseis trace, the seismic wavelet is the convolution of the zero-phase Klauder wavelet with the component minimum-phase wavelets. Thus the seismic wavelet occurring on a correlated vibroseis trace does not meet the minimum-phase requirement necessary for spiking deconvolution, and the final result of deconvolution is less than optimal. Over the years, this problem has been investigated and various methods of correction have been introduced. In essence, the existing methods of vibroseis deconvolution make use of a correction that converts (on the correlated trace) the Klauder wavelet into its minimum-phase counterpart. The seismic wavelet, which is the convolution of the minimum-phase counterpart with the component minimum-phase wavelets, is then removed by spiking deconvolution. This means that spiking deconvolution removes both the constructed minimum-phase Klauder counterpart and the component minimum-phase wavelets. Here, a new method is proposed: instead of being converted to minimum phase, the Klauder wavelet is removed directly. The spiking deconvolution can then proceed unimpeded as in the case of a dynamite record. These results also hold for gap predictive deconvolution because gap deconvolution is a special case of spiking deconvolution in which the deconvolved trace is smoothed by the front part of the minimum-phase wavelet that was removed.     

Spectral bandwidth enhancement of GPR profiling data using multiple-frequency compositing

The amplitude spectrum of ground penetrating radar (GPR) reflection data acquired with a particular antenna set is normally concentrated over a spectral bandwidth of a single octave, limiting the resolving power of the GPR wavelet. Where variously-sized GPR targets are located at numerous depths in the ground, it is often necessary to acquire several profiles of GPR data using antennas of different nominal frequencies. The most complete understanding of the subsurface is obtained when those frequency-limited radargrams are jointly interpreted, since each frequency yields a particular response to subsurface reflectivity. The application of deconvolution to GPR data could improve image quality, but is often hindered by limited spectral bandwidth.We present multiple-frequency compositing as a means of combining data from several frequency-limited datasets and improving the spectral bandwidth of the GPR profile. A multiple-frequency composite is built by summing together a number of spatially-coincident radargrams, each acquired with antennae of different centre frequency. The goal of the compositing process is therefore to produce a composite radargram with balanced contributions from frequency-limited radargrams and obtain a composite wavelet that has properties approximating a delta function (i.e. short in duration and having a broad, uniform spectral bandwidth).A synthetic investigation of the compositing process was performed using Berlage wavelets as proxies for GPR source pulses. This investigation suggests that a balanced, broad bandwidth, effective source pulse is obtained by a compositing process that equalises the spectral maxima of frequency-limited wavelets prior to summation into the composite. The compositing of real GPR data was examined using a set of 225, 450 and 900 MHz GPR common offset profiles acquired at a site on the Waterloo Moraine in Ontario, Canada. The most successful compositing strategy involved derivation of scaling factors from a time-variant least squares analysis of the amplitude spectra of each frequency-limited dataset. Contributions to the composite from each nominal acquisition frequency are clear, and the trace averaged amplitude spectrum of the corresponding composite is broadened uniformly over a bandwidth approaching two-octaves. Improvements to wavelet resolution are clear when a composite radargram is treated with a spiking deconvolution algorithm. Such improvement suggests that multiple-frequency compositing is a useful imaging tool, and a promising foundation for improving deconvolution of GPR data.     

基于矢量预测的SIMO系统混合相位地震子波提取方法研究

针对高阶统计量混合相位地震子波提取方法的局限性, 提出一种基于矢量预测的单输入多输出系统(SIMO)的混合相位地震子波提取方法. 该方法利用二阶循环平稳统计量包含的系统相位信息, 将CMP道集视为一个单输入多输出系统的输出, 利用道集中相邻两道或多道数据通过矢量预测来构建反子波计算的方程式, 进一步进行混合相位子波提取. 利用提取出的子波相位信息对CMP道集进行纯相位滤波, 能够替代相位校正技术, 并且利用获取的反褶积算子对CMP道集进行反褶积处理, 使提高分辨率后的不同道子波振幅、 频率、 波形相一致, 提高叠加的质量. 模型试算和实际资料处理结果表明, 文中方法适用于任意相位的子波提取及反褶积处理, 且处理精度较高, 具有较高的实际应用价值.      

Estimation of the seismic wavelet through homomorphic deconvolution and well log data: application on well-to-seismic tie procedure

Wavelet estimation and well-tie procedures are important tasks in seismic processing and interpretation. Deconvolutional statistical methods to estimate the proper wavelet, in general, are based on the assumptions of the classical convolutional model, which implies a random process reflectivity and a minimum-phase wavelet. The homomorphic deconvolution, however, does not take these premises into account. In this work, we propose an approach to estimate the seismic wavelet using the advantages of the homomorphic deconvolution and the deterministic estimation of the wavelet, which uses both seismic and well log data. The feasibility of this approach is verified on well-to-seismic tie from a real data set from Viking Graben Field, North Sea, Norway. The results show that the wavelet estimated through this methodology produced a higher quality well tie when compared to methods of estimation of the wavelet that consider the classical assumptions of the convolutional model.     

基于保角变换的双谱域相位谱估计方法及其在地震子波估计中的应用(英文)

地震子波估计是地震资料处理与解释中的重要环节,它的准确与否直接关系到反褶积及反演等结果的好坏。高阶谱(双谱和三谱)地震子波估计方法是一类重要的、新兴的子波估计方法,然而基于高阶谱的地震子波估计往往因为高阶相位谱卷绕的原因,导致子波相位谱求解产生偏差,进而影响了混合相位子波估计的效果。针对这一问题,本文在双谱域提出了一种基于保角变换的相位谱求解方法。通过缩小傅里叶相位谱的取值范围,有效避免了双谱相位发生卷绕的情况,从而消除了原相位谱估计中双谱相位卷绕的影响。该方法与最小二乘法相位谱估计相结合,构成了基于保角变换的最小二乘地震子波相位谱估计方法,并与最小二乘地震子波振幅谱估计方法一起,应用到了地震资料混合相位子波估计中。理论模型和实际资料验证了该方法的有效性。同时本文将双谱域地震子波相位谱估计中保角变换的思想推广到三谱域地震子波相位谱估计中。     

用遗传算法实现地震信号反褶积

遗传算法作为寻优手段具有全局优化和很好的稳定性．本文将遗传算法用于地震信号反褶积处理，与已往方法相比它具有更好的分辨率和稳定性我们采用Bernoulli－Gaussian模型和ARMA模型分别描述地震反射系数序列和地震子波，用最大似然和最小预测误差准则分别构造用于估计反射系数序列和地震子波的目标函数，用遗传算法优化目标函数，以实现地震信号反褶积．     

SIGNAL ADJUSTMENT OF VIBROSEIS AND IMPULSIVE SOURCE DATA*

In certain areas continuous Vibroseis profiling is not possible due to varying terrain conditions. Impulsive sources can be used to maintain continuous coverage. While this technique keeps the coverage at the desired level, for the processing of the actual data there is the problem of using different sources resulting in different source wavelets. In addition, the effect of the free surface is different for these two energy sources. The approach to these problems consists of a minimum-phase transformation of the two-sided Vibroseis data by removal of the anticipation component of the autocorrelation of the filtered sweep and a minimum-phase transformation of the impulsive source data by replacement of the recording filter operator with its minimum-phase correspondent. Therefore, after this transformation, both datasets show causal wavelets and a conventional deconvolution (spike or predictive) may be used. After stacking, a zero-phase transformation can be performed resulting in traces well suited for computing pseudo-acoustic impedance logs or for application of complex seismic trace analysis. The solution is also applicable to pure Vibroseis data, thereby eliminating the need for a special Vibroseis deconvolution. The processing steps described above are demonstrated on synthetic and actual data. The transformation operators used are two-sided recursive (TSR) shaping filters. After application of the above adjustment procedure, remaining signal distortions can be removed by modifying only the phase spectrum or both the amplitude and phase spectra. It can be shown that an arbitrary distortion defined in the frequency domain, i.e., a distortion of the amplitude and phase spectrum, is noticeable in the time section as a two-sided signal.     

Estimation of an optimal mixed-phase inverse filter

Inverse filtering is applied to seismic data to remove the effect of the wavelet and to obtain an estimate of the reflectivity series. In many cases the wavelet is not known, and only an estimate of its autocorrelation function (ACF) can be computed. Solving the Yule-Walker equations gives the inverse filter which corresponds to a minimum-delay wavelet. When the wavelet is mixed delay, this inverse filter produces a poor result.
By solving the extended Yule-Walker equations with the ACF of lag α on the main diagonal of the filter equations, it is possible to decompose the inverse filter into a finite-length filter convolved with an infinite-length filter. In a previous paper we proposed a mixed-delay inverse filter where the finite-length filter is maximum delay and the infinite-length filter is minimum delay.
Here, we refine this technique by analysing the roots of the Z -transform polynomial of the finite-length filter. By varying the number of roots which are placed inside the unit circle of the mixed-delay inverse filter, at most 2 ^α different filters are obtained. Applying each filter to a small data set (say a CMP gather), we choose the optimal filter to be the one for which the output has the largest L ^p -norm, with p =5. This is done for increasing values of α to obtain a final optimal filter. From this optimal filter it is easy to construct the inverse wavelet which may be used as an estimate of the seismic wavelet.
The new procedure has been applied to a synthetic wavelet and to an airgun wavelet to test its performance, and also to verify that the reconstructed wavelet is close to the original wavelet. The algorithm has also been applied to prestack marine seismic data, resulting in an improved stacked section compared with the one obtained by using a minimum-delay filter.     

An Approximation to Sparse-Spike Reflectivity Using the Gold Deconvolution Method

Wiener deconvolution is generally used to improve resolution of the seismic sections, although it has several important assumptions. I propose a new method named Gold deconvolution to obtain Earth’s sparse-spike reflectivity series. The method uses a recursive approach and requires the source waveform to be known, which is termed as Deterministic Gold deconvolution. In the case of the unknown wavelet, it is estimated from seismic data and the process is then termed as Statistical Gold deconvolution. In addition to the minimum phase, Gold deconvolution method also works for zero and mixed phase wavelets even on the noisy seismic data. The proposed method makes no assumption on the phase of the input wavelet, however, it needs the following assumptions to produce satisfactory results: (1) source waveform is known, if not, it should be estimated from seismic data, (2) source wavelet is stationary at least within a specified time gate, (3) input seismic data is zero offset and does not contain multiples, and (4) Earth consists of sparse spike reflectivity series. When applied in small time and space windows, the Gold deconvolution algorithm overcomes nonstationarity of the input wavelet. The algorithm uses several thousands of iterations, and generally a higher number of iterations produces better results. Since the wavelet is extracted from the seismogram itself for the Statistical Gold deconvolution case, the Gold deconvolution algorithm should be applied via constant-length windows both in time and space directions to overcome the nonstationarity of the wavelet in the input seismograms. The method can be extended into a two-dimensional case to obtain time-and-space dependent reflectivity, although I use one-dimensional Gold deconvolution in a trace-by-trace basis. The method is effective in areas where small-scale bright spots exist and it can also be used to locate thin reservoirs. Since the method produces better results for the Deterministic Gold deconvolution case, it can be used for the deterministic deconvolution of the data sets with known source waveforms such as land Vibroseis records and marine CHIRP systems.     

Lp-NORM DECONVOLUTION1

The purpose of deconvolution is to retrieve the reflectivity from seismic data. To do this requires an estimate of the seismic wavelet, which in some techniques is estimated simultaneously with the reflectivity, and in others is assumed known. The most popular deconvolution technique is inverse filtering. It has the property that the deconvolved reflectivity is band-limited. Band-limitation implies that reflectors are not sharply resolved, which can lead to serious interpretation problems in detailed delineation. To overcome the adverse effects of band-limitation, various alternatives for inverse filtering have been proposed. One class of alternatives is L_p-norm deconvolution, L₁norm deconvolution being the best-known of this class. We show that for an exact convolutional forward model and statistically independent reflectivity and additive noise, the maximum likelihood estimate of the reflectivity can be obtained by L_p-norm deconvolution for a range of multivariate probability density functions of the reflectivity and the noise. The L_∞-norm corresponds to a uniform distribution, the L₂-norm to a Gaussian distribution, the L₁-norm to an exponential distribution and the L₀-norm to a variable that is sparsely distributed. For instance, if we assume sparse and spiky reflectivity and Gaussian noise with zero mean, the L_p-norm deconvolution problem is solved best by minimizing the L₀-norm of the reflectivity and the L₂-norm of the noise. However, the L₀-norm is difficult to implement in an algorithm. From a practical point of view, the frequency-domain mixed-norm method that minimizes the L₁norm of the reflectivity and the L₂-norm of the noise is the best alternative. L_p-norm deconvolution can be stated in both time and frequency-domain. We show that both approaches are only equivalent for the case when the noise is minimized with the L₂-norm. Finally, some L_p-norm deconvolution methods are compared on synthetic and field data. For the practical examples, the wide range of possible L_p-norm deconvolution methods is narrowed down to three methods with p= 1 and/or 2. Given the assumptions of sparsely distributed reflectivity and Gaussian noise, we conclude that the mixed L₁norm (reflectivity) L₂-norm (noise) performs best. However, the problems inherent to single-trace deconvolution techniques, for example the problem of generating spurious events, remain. For practical application, a greater problem is that only the main, well-separated events are properly resolved.     

L1 norm inversion method for deconvolution in attenuating media

In order to perform a good pulse compression, the conventional spike deconvolution method requires that the wavelet is stationary. However, this requirement is never reached since the seismic wave always suffers high‐frequency attenuation and dispersion as it propagates in real materials. Due to this issue, the data need to pass through some kind of inverse‐Q filter. Most methods attempt to correct the attenuation effect by applying greater gains for high‐frequency components of the signal. The problem with this procedure is that it generally boosts high‐frequency noise. In order to deal with this problem, we present a new inversion method designed to estimate the reflectivity function in attenuating media. The key feature of the proposed method is the use of the least absolute error (L1 norm) to define both the data and model error in the objective functional. The L1 norm is more immune to noise when compared to the usual L2 one, especially when the data are contaminated by discrepant sample values. It also favours sparse reflectivity when used to define the model error in regularization of the inverse problem and also increases the resolution, since an efficient pulse compression is attained. Tests on synthetic and real data demonstrate the efficacy of the method in raising the resolution of the seismic signal without boosting its noise component.     

SEISMIC PHASE UNWRAPPING: METHODS,RESULTS, PROBLEMS1

Six known methods of seismic phase unwrapping (or phase restoration) are compared. All the methods tested unwrap the phase satisfactorily if the initial function is a simple theoretical wavelet. None of the methods restore the phase of a synthetic trace exactly. An initial validity test of the phase-unwrapping method is that the sum of the restored wavelet phase spectrum and the restored pulse-trace phase spectrum (assuming the convolutional model for the seismic trace) must be equal to the restored phase spectrum of the synthetic trace. Results show that none of the tested methods satisfy this test. Quantitative estimation of the phase-unwrapping accuracy by correlation analysis of the phase deconvolution results separated these methods, according to their efficiency, into three groups. The first group consists of methods using a priori wavelet information. These methods make the wavelet phase estimation more effective than the minimum-phase approach, if the wavelet is non-minimum-phase. The second group consists of methods using the phase increment Δø(Δω) between two adjacent frequencies. These methods help to decrease the time shift of the initial synthetic trace relative to the model of the medium. At the same time they degrade the trace correlation with the medium model. The third group consists of methods using an integration of the phase derivative. These methods do not lead to any improvement of the initial seismic trace. The main problem in the phase unwrapping of a seismic trace is the random character of the pulse trace. For this reason methods based on an analysis of the value of Δø(Δω) only, or using an adaptive approach (i.e. as Δω decreases) are not effective. In addition, methods based on integration of the phase derivative are unreliable, due to errors in numerical integration and differentiation.