地震子波提取方法综述

2.1.5 ARMA最小相位系统辨识 如果将地震数据表示成如下 ARMA 过程:
3
p
q
∑ ∑ x(n) + ai x(n − i) = e(n) + bje(n − j) (8)
i =1
j =1
其中 x(n) 为地震数据， e(n) 为地层反射系
数， ai 、 bj 为系统参数，可以得到：
∞
∑ Rx (τ ) = σ 2 h(i)h(i +τ )
(4)
其中 rbb (τ ) 为子波的自相关， rxx (τ ) 为地
震道的自相关，对其做 Z 变换：
M
∑ Rbb (z) = B(z)B(z−1) = Rxx (z) = rxx (τ ) ⋅ z−τ (5) τ =−M
对上式因式分解，一般可以得到 2M 个互为倒 数的零点，设单位圆内的零点为
zm ,
1.2 谱除法
这种方法遇到的一个问题是可能被零除， 因为反射系数谱中有零值或接近于零的值，习 惯上是将反射系数序列谱白噪化，以避免这种 现象的发生。
1.3 线性反演方法
通过正演模型计算出合成地震记录，然后 与实际地震记录比较，逐步修正子波，直到合 成记录与实际观测记录完全匹配为止。方法最 小化如下目标函数：
地震子波提取方法综述
杨培杰，印兴耀
中国石油大学（华东）地球资源与信息学院, 山东 东营 257061
摘要:如何准确提取地震子波是储层地球物理反演中一个非常关键的问题。地震子波提取方法主 要分为确定性子波提取和统计性子波提取两种，确定性子波提取方法首先利用测井资料计算出 反射系数序列，然后结合井旁地震道由褶积模型求出地震子波；统计性子波提取方法是通过地 震道自身来估计子波，它不需要测井资料的参与。文章对这两种方法进行了详细的总结，并指 出了它们的优点与不足；同时提到了一种新的可用于准确提取地震子波的子空间分解法 关键词：确定性子波提取，统计性子波的提取，混合相位子波，高阶统计量，子空间法
A(ω) ，待求的相位谱为ϕ(ω) ，则：
ϕ(ω) = −HT[A(ω)] = − ln A(ω) ∗ 1 πω
(3)
最小相位谱也可以用 Wold-Kolmogorov 公
式得到，然后与子波振幅谱组合再经过反傅立
叶变换得到最小相位子波。
2.1.2 Z 变换法
在反射系数为白噪声的假设条件下：
rbb (τ ) = rxx (τ )
(9)
i=0
线性系统的冲击响应就是冲击信号激励
该系统时的输出响应，故由 ARMA 过程的定
义式:
p
q
∑ ∑ aih(n − i) = bkδ (n − k) = bn (10)
i=0
k =0
综合上面两个公式可以得到下面的方程
p
∑ Rx (l) + ai Rx (l − i) = 0, ∀l > q (11) i =1
频率的关系为：Pr ( f ) = cf −β ，经过整理，我
们可以得到下面的公式：
Pw ( f ) = bPx ( f ) f β (7)
对上式两边开方并取绝对值得到：
|W ( f ) |=| bPx ( f ) f β | ，然后应用同态滤波的方 法结合希尔伯特变换得到相位谱，进而求取地 震子波。
2.1 最小相位子波提取[3][10]
2.1.1 自相关法
在反射系数为白噪声的假设条件下，地震 记录的自相关等于子波自相关，通过地震记录 的自相关得到子波的振幅谱，对于子波是零相 位的情况下，可以直接由地震记录的振幅谱得 到子波的频谱，由傅立叶反变换得到子波的时 间域形式。当子波为最小相位时，可以由希尔 伯特变换法得到相位谱，记子波的振幅谱为
通过地震道的四阶矩，应用脉冲反褶积的 方法求得一个反子波，然后由反子波得到最小 相位子波得振幅谱和相位谱，最终得到最小相 位子波，该方法可以有效的抑制地震道中的高 斯白噪声和高斯色噪声。
2.1.4 分形地震子波提取[11][12]
在假设地震记录所反映的波阻抗函数具 有分形特性条件下，得到反射系数的功率谱与
统计性子波提取方法是通过地震道自身 来估计子波，又可分为基于二阶统计量和基于 高阶统计量两种方法。二阶统计量方法首先由
1
Robinson(1975)提出，它基于这样的假设，即 地震子波是时不变的，地下的反射时具有白噪 谱的随机序列，则观测到的地震道的自相关就 给定了地震子波的自相关的一个估计，也就是 已知了地震子波的振幅谱，对于子波的相位 谱，则必须给出一定的假设，如假设地震子波 是零相位、最小相位、最大相位，而实际上地 震子波是一种混合相位的，因此，基于二阶统 计量的自相关统计的方法提取的子波也是不 准确的。80年代后期以来，许多学者开始使用 高阶统计方法来解决地震子波估计的问题，这 些方法大多源于60年代发展起来的累积量和 多谱理论[2]，T. Matsuoka和T.J. Ulrych(1984) 最早将它们用于混合相位地震子波的估计，高 阶统计量地震子波提取的新思想是由 G.D.Lazear(1993)和D.R.Velis(1996)提出的，他 们将非高斯信号处理中的四阶累积量用于子 波的估计，为解决混合相位子波估计问题提供 了一条全新的思路。
(m
=
1,L,
k)
，及
1 zm
,
(m
=
k
+ 1,L ,
M
)
，
则：
∏ ∏ M
k m=1
(
z−zm
M
) (z−
m=k+1
1 zm
)
B(Z) = zM
(6)
= b0 +b1z−1 +b2z−2 +L+bM z−M
由此可以得到最小相位地震子波
b(t) = (b0 ,b1,L, bM ) 。
2.1.3 四阶矩法
Subspace method
前言
1地震子波是地震勘探中一个非常关键的问题 [1]，在正演问题中，需要通过波动方程或褶积 模型结合地震子波来形成正演模拟地震数据， 在反演和反褶积问题中，也需要通过地震道提 取一个的子波，不同的子波往往对反演结果会 有不同的影响。地震子波是一段具有确定的起 始时间、能量有限且有一定延续长度的信号， 它是地震记录中的基本单元，一般认为，地震 震源激发时所产生的地震波仅是一个延续时 间极短的尖脉冲，随着尖脉冲在粘弹性介质中 传播，尖脉冲的高频成分很快衰减，波形随之
除此之外，还又很多其它的确定性地震子 波提取方法，这里就不再一一赘述。
2 统计性地震子波提取
2
从盲信号处理的角度出发，统计性子波提
取应该属于盲系统辨识[8][9]的范畴，如图 1 所
示，盲系统辨识(Blind System Identification，
BSI)就是指在输入信号R和系统函数W都未知
Survey on seismic wavelet extraction
YANG Pei-jie, YIN Xing-yao
(College of Earth Resources and Information, China University of Petroleum (East China), Dongying 257061，Shandong Province，China)
假定要估计的地震子波为θ ，为了进行贝 叶斯推断，首先要来设定θ 的先验分布 p(θ ) ，
然后将得到的观测信息融入先验信息wk.baidu.com改进
先验分布，假定地震数据x(它是与θ 有关的)， 用 p(x |θ ) 表示其条件密度函数。有了先验分
布和条件分布，就可以通过贝叶斯公式得到参
数θ 的后验分布：
p(θ | x) = h(x,θ ) ∝ p(θ ) p(x |θ ) (2) p(x)
minf (w) = (r*w-d )(r*w-d )+Q(w) (1) 其中，w 为待求的地震子波，r 为地层反 射系数，d 为地震数据，Q(w) 为对子波的先验 约束。
1.4 贝叶斯法[6]
该方法运用贝叶斯公式将地震子波、反射 系数、地震道噪声的先验信息与似然函数相结 合从而得到各个参数的后验概率密度。
的情况下，仅仅根据系统的输出D来辨识系统
函数W。
输入 R
输出 D
系统函数 W
图1 输入和系统响应未知的盲系统辨识示意图
统计性地震子波提取与盲系统辨识问题 类似，如果将地层反射系数看作输入 R，将系 地震子波看作是系统函数 W，统计性地震子波 提取就是在地层反射系数和地震子波都未知 的情况下，仅仅根据观测到的地震记录 D 来估 计地震子波。
Abstract: How to extract seismic wavelet accurately is a very key problem while seismic inversion. Wavelet extraction method are divided into deterministic wavelet extraction and statistical wavelet extraction mainly, deterministic wavelet extraction utilizes well logging information to calculate out the reflection coefficient at first, then combine with seismic data near the well to extract wavelet by using convolution model; statistical wavelet extraction only utilizes seismic data alone to extract seismic wavelet, it does not need participation of the logging information. This article carried on detailed summary to these two kinds of methods, and pointed out their advantages and disadvantages; meanwhile, a new wavelet extraction method named subspace method is mentioned Key words: Deterministic wavelet extraction, Statistical wavelet extraction, Mix-phase wavelet, High-order statistics,
这一法方程就是著名的修正 Yule-Walker
方程，特别地，对于严格 AR(p)过程，上式简
化为：
p
∑ Rx (l) + ai Rx (l − i) = 0, i =1
∀l > 0 (12)
通过求解以上方程组，就可以确定 AR 参
数 a1,L, ap 。
1 确定性地震子波提取
确 定 性 子 波 提 取 方 法 [3][4][5] 是 利 用 声 波 测 井和密度测井资料，首先计算出反射系数序 列，然后结合井旁地震道由褶积模型求出地震 子波：
1.1 维纳滤波法
它是在时间域内通过解一线性方程组估计 子波，实际计算过程中只能截取一段地震记录 来估计子波，这就会产生截断误差。
通过最大后验估计(maximum posteriori, MAP)或马尔可夫链-蒙特卡罗(Markov chain Monte Carlo，MCMC)模拟方法得到地震子波。
1.5 循环迭代法[7]
方法考虑到地震子波在空间变化的特点， 首先由多道相关法提取初始的地震子波的振 幅谱，然后结合测井资料确定初始地震子波的 相位谱，然后根据离散反演理论迭代求取精细 的井旁地震子波。
基金项目：中国石油大学优秀博士学位论文培育资助 项目(B2007-02)
增长，便形成了地震子波，一个地震子波一般 有二至三个相位的延续长度，大约有90ms左 右，然后以地震子波的形式在地下传播。
地震子波提取的基本框架是褶积模型，也 就是包含子波以及反射系数序列再加上噪声 的地震道，提取方法主要包括两大类，第一类 是确定性子波提取方法(deterministic wavelet extraction methods)，第二类是统计性子波的提 取方法(statistical wavelet extraction methods)， 确定性子波提取方法指的是利用测井资料首 先计算出反射系数序列，然后结合井旁地震道 由褶积模型求出地震子波。