1、第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,反褶积又称反滤波。为了消除大地滤波及接收系统滤波对地震数据的影响而作出的滤波处理。反滤波本质上是一种频率滤波。从数学上看，它是一种褶积运算，故称反褶积。,3.1、反褶积及褶积模型,一、反褶积的概念,反褶积处理：是常用处理方法之一。可以用于叠前和叠后，也可以多次使用。作用：压缩地震子波，提高分辨率。 可以压制多次波和短周期鸣震等干扰，提高地震资料信噪比。,震源爆炸使地下介质形成三个区域：,地层
2、对震源脉冲的改造作用，相当于对它进行了一次低通滤波，此滤波器常称为大地滤波器。,震源爆炸产生尖脉冲传播到弹性区起始边界时，已经变成了有一定延续时间的稳定波形地震子波。,假设震源脉冲在地下介质中传播未受大地改造，脉冲信号入射到分界面、反射信号返回地面，被检波器接收、传输到仪器被记录下来。如果接收系统未对震源脉冲进行改造，则地震记录为反射系数序列：,反褶积就是要获得未经系统作用的地震波形。,地震子波 同震源子波 ，其概念是有区别的，它与许多因素有关。根据地震波传播过程中影响因素的不同，地震子波可描述为：,式中,而,干扰波是由非激发干扰（次生） 、背景噪声及规则（或称相干）（由激发产生）干扰 叠加而
3、成:,规则干扰 分两类：一类与地质构造有关，包括多次波、转换波、绕射波、伴随波、折射波、瑞利波、勒夫波和斯通利波等,这类波在特定的条件下可转化为有效波;另一类与地质构造无关,如水中震鸣、气泡效应、地表及海面散射等 (也包括地下震鸣、薄层微曲多次波)。实际处理时,要根据不同的勘探情况,分别对待。,反褶积的关键是如何设计一个反滤波器去抵消另一个滤波器的作用。设计反滤波器的方法：,由已知地震子波计算反褶积算子，称确定性反褶积，主要用于去除记录系统的响应、海上震源子波反褶积等方面；通过统计方法求取最佳反褶积算子，如脉冲反褶积、预测反褶积等。,二、褶积模型,理想模型：,加噪模型：,其中：,（3-1）,（
4、3-2）,反褶积的假设条件： （1）地下地层是水平层状介质； （2）地震波是垂直入射反射的平面波； （3）地震子波在传播过程中保持波形不变； （4）地震记录中无噪声； （5）地震子波已知； （6）反射系数序列为白噪序列； （7）地震子波是最小相位的。若假设条件与实际不吻合，势必会造成褶积模型与实际地震记录存在一定差异。,改进模型：,海上“特征反褶积”模型,其中：,沙漠地区可控震源地震记录模型,其中：,仿真褶积模型,其中：,为了把地震子波压缩成尖脉冲（必需去掉大地滤波器的作用），使地震记录变为反射系数序列，出现了各种反褶积方法，而实际处理结果往往不如人愿。其原因有三：地震记录已知，地震子波未知，
5、求反射系数序列，必须有若干假设条件限定解的唯一性，否者是多解的；假设条件与实际情况越接近，反褶积效果越好。,反褶积方法依赖地震记录的褶积模型，模型中地震子波是大地滤波器的脉冲响应，而大地滤波的作用复杂，模型不太可靠。只有先彻底解决正演问题，才能使反褶积得到发展。反褶积方法可能会提高噪声水平，有必要同时发展提高分辨率及信噪比的方法。反褶积方法很多，有些（如最大熵、卡尔曼、时变Q等）未能在常规处理中获得一席之位。,反射系数剖面,地震剖面,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3
6、.7、同态反褶积 3.8、地表一致性反褶积,3.2 反滤波 一、反滤波的概念 1、概念,假定地震记录不含干扰，即,2、反子波,对应的频域形式,则可得到,写成时域形式为：,反子波与子波褶积为：,由子波和反射系数求地震记录，是一褶积过程（正演）；已知反子波和地震记录求反射系数，称为反褶积或反滤波。,二、地震子波的求取,确定性反褶积，需已知子波。故先讨论子波求取方法，有5种方法： 直接观测法（适用于海上）；自相关法；多项式求根法；利用测井资料求子波；对数分解法。确定性反褶积处理步骤：先提子波，再求反子波，然后进行反褶积。,2、自相关法,选择一段质量较高的地震记录，时窗长度为T：,其Z变换为,假设反射
7、系数是白噪声序列，其z变换为 则 的自相关 的z变换：,地震记录的z变换为,地震记录自相关 的z变换为,将 代入，有：,由于 都是实数序列，所以有：,因此有：,也有：,未知，现在来确定它,（3-19）,（3-18）,假如地震子波是最小相位的物理可实现序列，则其z变换为：,由物理可实现性知：当 时， , 对下式,（3-19）,两端取对数，有：,令,因而得到,（3-20）,根据复变函数理论,其中C为常量，和,因而得到：,因为 是实数，由（3-20）式知 ，于是可得C,（3-20）,对上求希尔伯特变换，即可求出相位谱：,令,（3-19）,求出相位谱后，对下式,求付立叶反变换，得最小相位子波,（3-2
8、3）,复杂内容简单化：,假如地震子波是零相位的，下式,（3-19）,中的相位谱为零，即 ，因此有,零相位地震子波：,（3-32）,3、多项式法,选择一段质量较高的地震记录,设反射系数为白噪声序列，则记录自相关与子波自相关等价，即,褶积模型：,令 ，则,（3-34）,将上式两端乘以 ，则有：,由于,显然， 应有2M个根。鉴于系数均为实数，所以2M个根是M对互为倒数的，即若则另一根为：,根据这M对根在单位圆内、外的位置，可以组成2M个不同相位的地震子波，其中必有一个是最小相位，一是最大相位的。根据“最小相位序列z域零点在单位圆内”这一特点，选出模小于1的根，便可组成最小相位子波，其z变换为：,令z
9、=0,得,由此得最小相位子波,例如，已知记录的自相关,其z变换为：,两端乘以 z2 得：,求 的根，有,选出模小于1的根，求,得最小相位子波为：,4、用测井资料求取子波这种方法要求有较好的声波测井和密度测井资料，并在井旁有质量较高的地震记录,首先将声波时差 转化为声波速度 设声波单位为 ，则,然后，进行深时转换， 是双层旅行时,然后，计算反射系数对井旁地震记录和反射系数进行傅氏变换后，可得到子波的频谱对 傅氏反变换就得到地震子波,即：,5、对数分解法,时域模型：,此法不需假设反射系数是白噪声，不需假设子波是最小相位。,频域模型：,对频域模型两端取对数，则将子波与反射系数分离开来,称 为对数谱。
10、,用付立叶反变换对数谱的时间信号：,由于 分布在时间轴原点附近； 分布离原点较远区域。若二者分离较好，则可用低通滤波将 分离出来，便可求出子波 。,即先付氏正变换，再取指数，然后进行付氏反变换,由于很难确定对数谱 在时间轴上的分布区域，故用N道地震记录来求取。,求平均值：,由于 ，所以,因此可获得较准确的子波。,三、反滤波的实现,求出地震子波 后，可用付氏变换或z变换或最小平方法求反子波 ，然后对地震记录 做褶积运算得 。,如,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、
11、同态反褶积 3.8、地表一致性反褶积,3.3 最佳维纳滤波及最小平方反褶积,一、最佳维纳滤波,维纳滤波即最小平方滤波, 它是使实际输出与期望输出，在误差能量最小条件下，求滤波因子的方法，所代表滤波器是最佳滤波器。,求解关系输入信号滤波因子实际输出期望输出误差能量,最小平方最小,滤波方程,写成矩阵形式,当反射系数为白噪声时，记录和子波自相关等价，也就不用求子波了。,表3-1 常用的期望输出,根据滤波目的设定期望输出,二、最小平方反褶积,将最佳维纳滤波原理应用于反褶积问题，就是最小平方反褶积方法。地震记录为：,地震记录含随机噪声 后,褶积模型变为：,设计滤波因子为：,输出为：,误差的平方和为：,式
12、中：,表明：地震记录的自相关等于子波自相关与噪音自相关的和,表明：地震记录与期望输出的互相关等于期望脉冲因子与子波的自相关,写成矩阵形式,（3-76）,（3-75）,地震记录含噪声后,求得的反子波与不含噪时会有差别(托布里兹矩阵主对角线元素上加了噪声的自相关值e),但这一点并不影响反褶积的效果,反而可增加方程组求解的稳定性。有时为了实际需要,要人为加进一些噪声,这就是预白化问题。,例1:,（1）求实际输出序列,（输入：最小相位子波，期望输出：零延迟脉冲）,按原理求,（2）求滤波因子,例2:,（输入：最大相位子波，期望输出：任意延迟脉冲）,按方程求,三、最小平方反褶积中的预白化处理,（不讲）,对
13、（3-43）式描述的是频谱形式, 对其改造为（3-43）为了解决带限问题,在地震信号的功率谱P()中,从低频到高频统一加一白噪声。,2.预白化处理,（3-44）,当反射系数为白噪声时，子波的自相关=地震记录的自相关。,从以上方程可以看出：,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,3.4 脉冲反褶积,一、脉冲反褶积原理,当期望输出为 时， 除 外，其余全为零。,方程（矩阵形式）：,（3-80）,（3-79）,脉冲反褶积（对应3-80
14、式）要求输入子波为最小相位的，实际中往往是非最小相位的，因此要作某种优化处理。,最小平方反滤波，对子波相位已无要求，可以任意。,反滤波因子在理论上应为无穷项，实际只能取有限项，并且其主要部分,反滤波因子的形状和位置要根据地震子波的形状来决定。,详情见书p78,最小延迟，m0=0，m,最大延迟，m0=-n-m,混合延迟，m0=-m1,二、参数选择（效果好的标准是什么？子波是否被压缩）,（1）反滤波因子长度m，通过试验来选择；,（2）相关时窗长度m+n，常为2m记录长度；,（3）稳定常数，一般取 的百分数。干扰小，取0.0050.01 ；干扰大，取0.020.05 。,下图是延迟脉冲反褶积效果分析
15、图,（误差能量）,延迟时间,输入为非最小相位子波，期望输出,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,3.5 预测反褶积,一、预测滤波原理,预测问题是已知某个物理量的过去值和现在值,通过对已知信息加工处理来获得未来某个时刻的预测值。其数学描述如下:,设 为现在值， 为过去值, 如果定义预测步长为 ，用现在值和过去值来预测将来 时刻的预测值 ，即,（3-83）,使预测值与实际未来值的误差预测误差,最小，按最小平方原理，得如下方程：,（
16、3-86）,矩阵形式为：,（3-87）,二、预测反褶积原理,预测反褶积是要从含有多次波的地震记录中，根据多次波具有周期性特点，预测出多次波来，用前者减去后者得含一次波的地震记录（预测误差）。即：,地震记录=一次波（不可预测）+多次波（可预测）,地震记录-多次波（可预测） = 一次波（不可预测）,1、预测反褶积原理,设地震子波满足最小相位条件，反射系数为白噪声,褶积模型为,则 时刻的输出值为：,含未来时的信息,含过去、现在时的信息,比较如下两式：,可以看出: 是预测值 (即鸣震干扰), 是预测误差 (一次反射波).,设 即相当于取了子波的前部分,则式变为,式表明：将一个子波的前部分与反射系数的褶
17、积就得到了预测误差（一次反射波）。这种方法也压缩了子波的长度，从而提高地震资料的分辨率。,两步法实现过程：,设输入序列为5点序列： ，预测步长 ，期望输出：,、求预测因子,、求预测值,、求预测误差,即先做预测反滤波，再求预测误差，故称两步法。能否直接求预测误差？必须找出预测反滤波因子和预测误差反滤波因子的关系,预测反褶积数学描述：,其中 为预测误差反滤波因子，它为：,预测反褶积使子波长度变为 ，故也称其为子波切除反褶积。当 时，预测反褶积变为脉冲反褶积。,2、预测因子的求取,即：,输入信号预测因子实际输出期望输出,3、计算举例,已知地震记录存在虚反射,其时间序列为,预测反褶积来消除虚反射。,解
18、：设预测因子长 ，预测步长 。,(2)解方程,得预测因子：,(3) 求虚反射,即预测值:,(4)预测误差,消除虚反射后的那一段地震数据,(5)压制虚反射后的地震记录,周期为5,因周期为5，故预测步长取5,或,MATLAB 程序：,三、预测反褶积压缩反射脉冲（自学） 结论：,四、参数选择,（1）预测步长，关键参数，效果与此有关； （2）预测因子长度，要适宜，小，波形尾部波动；大，增加运算量，相邻反射及噪声影响； （3）预白化量，增加求解的稳定性，需根据噪声水平来定。后有三图，以助说明。,预测因子长度128ms，预测步长2ms，效果最好，其余有延续时间（非脉冲信号）。,（a）为脉冲序列，（b）为合
19、成记录,预测步长2ms,预测因子长度128ms，效果好。小长度因子，尾部波动明显,预白化量0.1，效果较好；预白化量增大，尾部波动增大，效果变差。,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,3.6 子波整形反褶积子波的相位与分辨率子波的相位通常有三种，即最小相位，混合相位和最大相位，这些子波是单边的物理可实现信号海上勘探和陆上勘探爆炸震源产生的地震子波都接近最小相位，所以，经常假设和讨论最小相位还有一种相位的子波是零相位子波，虽说是
20、一种非因果信号，是物理不可实现的，但在数字滤波、反褶积和反演中经常用到。,为什么要用到零相位子波，基于以下两个考虑：陆上可控震源子波因是通过自相关处理得到的，所以是零相位的；零相位子波，比其他相位子波的分辨率高,子波与反子波的时域分布的特点对于褶积公式假设子波为有限长度，从理论上来说，我们可以得到无限长的子波利用变换，我们可以推导出有限长子波与所求反子波的关系,.相位对反褶积精度的影响,在前面的讨论中，为了讨论问题方便，都假设了子波为最小相位实际上子波有不同的相位在反褶积中，其他条件相同，相位不同，反褶积的精度是不同的对于最小平方反褶积，滤波方程为,我们分两种相位的期望输出来讨论反褶积问题）期
21、望输出为（0，1）这时的方程为解得反子波为 ，反褶积输出为与期望输出（0，0，1）的误差能量为 ,从中我们也可以看到，此时的误差比最小相位是大多了 ）期望输出为（0,1,0） 这是的滤波方程为得到的反子波为 ，反褶积输出为这与期望输出(0,1,0)的误差能量 为 式中我们可以看到此时的反褶积比上面的精度要高,由此我们可以得到一个重要的结论：子波振幅相同时，最小相位子波对期望输出为零延迟的反褶积，误差最小；在子波为混合相位和最大相位时，期望输出的相位应与子波的相位相匹配，有一个最佳延迟，只有这样才能得到合适的反褶积结果,最小相位子波，期望输出波形为零延迟尖脉冲；最大相位子波，期望输出波形为最大延
22、迟尖脉冲；混合相位子波，期望输出波形为零延迟尖脉冲与 最大延迟尖脉冲之间的非零延迟尖脉冲。,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,第三章 反褶积,3.1、反褶积及褶积模型 3.2、反滤波 3.3、最佳维纳滤波与最小平方反褶积 3.4、脉冲反褶积 3.5、预测反褶积 3.6、子波整形反褶积 3.7、同态反褶积 3.8、地表一致性反褶积,3.8 地表一致性反褶积,目的：消除由于近地表变化对地震子波波形的影响。,地表一致性谱分解（Ta
23、ner 1981）认为任意道的频谱是对应的炮点、接收点、共中心点和炮检距的频谱响应的乘积，对其取对数，可分解炮点和接收点的振幅响应和相位响应。振幅响应被用于振幅一致性补偿，相位响应被用于时差校正。,地表一致性谱分解的基础是褶积模型和预测原理，并扩展为多道平均，把等效子波分解成多个分量。,地表一致性模型：,其中：,求 频谱 ：,其振幅谱和相位谱分别为：,对振幅谱取对数:,设 是实际地震道的振幅谱，根据地表一致性的假设 ，一个地震道由相应分量组成，而每一分量对相关的地震记录的贡献是一样的，则实际振幅谱与模型估计振幅谱之差的最小平方和为：,使误差能量达到最小，得到关于4个分量的方程组。,用迭代法求解
24、取每个分量，在最小相位条件下，估算它们的反算子，进行反褶积。,与单道反褶积的区别：炮点（或接收点）分量里仅含与观测点有关因素，对相关的所有道是相同的。,一方法原理褶积模型对于地震道有式中w(t)为一中和子波，其他符号意义同前考虑到地表四个因素：,共炮点,共接收点,共中心点,共炮检距点,（3-118）,式中,上式在频域式,分别写出它们的振幅谱和相位谱的形式,2.对数谱分析方法以上子波包括4个振幅谱和4个相位谱，假设w(t)为最小相位子波，这样只须讨论子波的振幅谱。为了计算方便，对（3-120）式两边取对数，得到于是，褶积关系变成相加关系，下标i代表第i炮，下标j代表第j道。简记对数振幅谱为 ，并
25、误差函数为（3-123）,这样在每个i、j上，可由（3-122）得到一个 ，它与实际地震道 有误差。令E最小，用Gauss-Seidel方法，可以求出4个对数振幅谱分量，再用反对数变换得到他们的振幅谱。反褶积因子的计算有了子波的振幅谱后，再假定子波为最小相位，可求出个反子波，对数据道相继采用这些反子波，就可以完成地表一次性反在褶积,二、地表一致性反褶积的应用反褶积之前的有关处理在反褶积之前，必须先进行去噪声处理，去除面波以及规则的干扰；同时要做好各种振幅恢复补偿处理时窗的选取为了得到各个频率分量的振幅谱和自相关函数，要控制剖面上地质构造的形态，划分时窗，每一道时窗按半个时窗互相重叠得到的各个分
26、时窗的自相关函数后，再做统计平均，作为托布里兹矩阵的自相关函数,反褶积因子的选取用以上个分量实现的反褶积，计算量十分巨大，是普通反褶积的倍在同一道集中，共中心位置分量来自相同地下界面，计算时可以忽略因此，通常情况下，只采用个分量对道进行反褶积但计算量还是很大的，实际使用时，不是对每一分量求反褶积因子，而是把个分量的振幅谱求和，求一个综合的反褶积因子另外，还可以在时间域，用迭代的方法求解,如果忽略检波点的因素，只考虑炮点和接收点的情况，地表一致性反褶积实际变成了两步法统计子波反褶积再退一步，如果只做炮点的反褶积，就变成一步反褶积，此时已不再称为一致性反褶积了效果分析图3-19是爆炸源的叠加地震剖
27、面分别做普通反褶积和地表一致性反褶积的结果比较可看出，地表一致性反褶积信噪比要比普通反褶积高，并且振幅变化均匀，特别是在.之间的浅层,地表一致性反褶积的主要目标是校正子波的振幅谱，所以，这种方法适合于地表条件变化大的地区但是真正方法并不是着重展宽频谱，所以，分辨率并不能提高很多,3.9 反Q滤波及谱白化,由于地层的吸收作用,地震波经地层传播后,能量被衰减损耗,频率变低,特别是深层,分辨率大大下降.因此,为恢复地震波原能量,必须做吸收补偿,即Q补偿或滤波.,1、Q因子的物理含义,波在介质中传播，产生的弹性能逐渐被介质吸收，最后转换为热能，其过程称为吸收。一个波长内，原地震波能量与传输所损耗的能量
28、之比为Q。,2、Q值的估算 (1)振幅包络法； (2)李氏经验公式：,3、反Q滤波的方法实现,设反Q滤波器为最小相位，振幅谱为,求出滤波因子，便可滤波。, Set Deconvolution Type to Minimum Phase Predictive Use an Operator Length (ms) of 250 Set the Prediction Distance (ms) to 96 Let White Noise (percent) be .1, the default Use one time gate with a Start Time (ms) of 700 and
29、an Interval(ms) of 400.,Predictive Deconvolution Processing,Parameter Options and Time Gate design for Deconvolution,Line 127 With 96 ms Predictive Deconvolution,Line 127 Without Post-Stack Processing Applied,Line 127 With 12 ms Predictive Deconvolution,Line 127 With 24 ms Predictive Deconvolution,P
30、ower Spectrum from line 127 with no decon.,Power spectrum of line 127 after minimum phase predictive deconvolution with a 36 ms gap,Power Spectrum from line 127 with no decon.,Power spectrum of line 127 after minimum phase predictive deconvolution with a 24 ms gap,Power Spectrum from line 127 with n
31、o decon.,Power spectrum of line 127 after predictive deconvolution with a 4 ms gap. This is an undesirable result!, Use .1% White Noise Use 1% White Noise Use 10% White Noise Use .1% White Noise with Zero Phase Spiking Deconvolution Type All other parameters remain the same.,Spiking Deconvolution Pr
32、ocessing,Line 127 Zero Phase Spiking Deconvolution.,Line 127 Minimum Phase Spiking Deconvolution.,Line 127 Power Spectrum Before Spiking Deconvolution,Line 127 Power Spectrum After Spiking Deconvolution, 10% Noise,Phase Correction (No Whitening-不白噪化), The Deconvolution Type is Phase Correction Only
33、Set Operator Length (ms) to 250 Use .1% White Noise Use one time gate with Start Time (ms) of 700 and Interval (ms)of 400,Line 127 After Phase Correction Deconvolution,Line 127 Before Phase Correction Deconvolution,Power spectrum of line 127 before phase correction deconvolution,Power spectrum of li
34、ne 127 after phase correction deconvolution.Note little change in frequency content,Spectral Balance谱均衡,Set the Spectral Balance Parameters: Number of Frequency Panels is 6 Scaling Time Window (250) Set the corner points （拐点）as follows:,Power spectrum of line 127 before spectral balancing,Power spectrum of line 127 after spectral balancing with 6 frequency panels,Line 127 with Phase Correction and Spectral Balance,
