反褶积学习基础课件

四川石油管理局地球物理勘探公司,欢迎各位光临！,欢迎各位光临！,地震资料处理中的反褶积处理,黄大云,2008年6月,地震资料处理中的反褶积处理黄大云,主要内容,1、在地震资料处理中为什么要做反褶积处理,2、预测反褶积的基本原理和计算方法,3、预测反褶积处理,4、地表一致性反褶积处理,5、子波整形反褶积处理,6、调谐反褶积,7、谱白化处理,8、反Q滤波,9、关于单道反褶积的讨论,主要内容,地震资料常规处理一般流程,数据输入,观测系统加载,坏炮坏道去除,噪声压制,振幅补偿,静校正,反褶积,叠加,剩余静校正,速度分析,动校正,偏移,反褶积-提高分辨率,地震资料常规处理一般流程数据输入观测系统加载坏炮坏道去除噪声,为什么要做反褶积处理?,为了说明在地震资料处理中要做做反褶积，有必要弄清褶积和反褶积的意义。,1、褶积的定义,褶积是一种数学运算的方式以及运算结果。定义如下：,两个函数x(t)和y(t)的褶积定义为：,在离散有限的情况下，积分变成以下求和形式：,我们通常用到的多为离散有限的情况。从以上公式可以看出，褶积就是先将其中一个函数（序列）反转过来再对应相乘并求和。即所谓的先褶后积，褶积的名称由此而来。,为什么要做反褶积处理?,2、地震记录的褶积模型,设震源发出的信号为 b(t),它遇到第一个到第n个反射界面的反射系数分别为g1、g2、gn,则检波器接收到各处的反射信号分别为g1.b(t-t1)、g2.b(t-t2)、gnb(t-tn)，地震记录x(t)为各反射信号之和，即：,上式表明：地震记录由地震信号和反射系数序列的褶积构成：,x(t)=g(t)*b(t),3、地震记录的分辨率,地震记录的分辨率由地震信号（地震子波）b(t)的延续长度和反射系数g(t)之间的距离决定。b(t)的延续长度越短，g(t)之间的距离越大，分辨率越高，反之分辨率越低。,2、地震记录的褶积模型,通常震源产生的信号（震源子波）是较短的，但它在传播过程中，因为大地的滤波作用会逐渐拉长，以至分辨率越来越低。为了提高分辨率，只有两种办法：加大g(t)之间的距离或者压缩b(t)的延续长度。,g(t)之间的距离是客观存在，显然我们无法也不应该去改变它。为了提高地震记录的分辨率，只有压缩地震子波b(t)的长度。理想的情况是将b(t)缩为单位脉冲函数。这时地震记录x(t)就是反射系数序列g(t)。把b(t)缩为单位脉冲函数的方法通常是用某种办法设计出一个算子，它与b(t)褶积的结果就是单位脉冲函数。由此我们得反褶积的定义。,通常震源产生的信号（震源子波）是较短的，但它在传播过程,4、反褶积的一般定义,反褶积就是去掉地震记录中大地的滤波作用的一种处理方法，所以反褶积也叫反滤波。它用的运算方法归根到底仍然是褶积。,但现在的反褶积已不局限于去除大地的滤波作用，凡是对地震子波进行改造的处理都叫它反褶积。,5、反褶积处理的目的,提高地震记录的分辨率是反褶积处理的目的之一，但对叠前反褶积而言，它却不是主要目的。叠前反褶积的主要目的是使地震子波波形一致，以便获得好的叠加效果。,4、反褶积的一般定义,未做反褶积,地表一致性反褶积后,未做反褶积地表一致性反褶积后,预测反褶积的原理和计算方法,主要内容,什么叫预测,预测的条件,预测滤波,预测反褶积,预测反褶积的计算,Omega系统的预测反褶积模块及其主要参数,因为Omega系统的几个重要的反褶积模块，如预测反褶积、地表一致性反褶积和调谐反褶积都用的是预测算法，所以有必要先说一说预测反褶积的原理和计算方法。,预测反褶积的原理和计算方法主要内容 因为Omega,什么叫预测?,预测就是根据过去和现在已发生的事实判定将来会出现的情况。,在数学上，对一个时间函数的预测是指该函数某一点的值用其前面若干个值的线性组合表示出来。这种预测称为线性预测。,什么叫预测? 预测就是根据过去和现在已发生的事实判定将来,预测的条件,并非所有事物都可线性预测。函数x(t)可线性预测的条件是:,x(t)为平稳随机过程，即它的统计特征：数学期望,和方差,是与时间无关的量，且自相关函数,r,xx,()只与时差有关 。,我们认为地震记录满足以上条件，因而可做预测。,预测的条件,预测滤波,在地震勘探中，我们认为地震记录是平稳随机过程，因而可以预测。,根据地震记录褶积模型的假设，地震记录x(t)由地震子波b(t)和地层反射系数g(t)的褶积构成：,我们先假定b(t)为一物理可实现的最小相位信号， g(t)为白噪序列。在时刻(t+),地震记录的振幅值可表示为：,在右端第二项中，令j=s-,上式变为：,记,设b(t)的反信号为a(t)，有a(t)*x(t)=a(t)*b(t)*g(t)=(t)*g(t)=g(t),因为b(t)为一物理可实现的最小相位信号，因此有：当t0),预测输出：,预测误差：,误差总能量：,预测反褶积的计算 预测反褶积计算的关键是求得预测滤波,选取c(s)，使Q达到最小。为此令,或,记,于是有：,选取c(s)，使Q达到最小。为此令,将以上方程写成矩阵形式就是：,以上方程的系数矩阵和左端的向量均由 x(t) 的自相关函数构成。该方程叫做预测方程，求解此方程，即得到最小平方意义下的预测滤波因子c(s)，用c(s) 对x(t) 滤波，若输出x(t+),就是预测滤波，若输出e(t+),就是预测反滤波或预测反褶积。,因为,所以预测反褶积算子为 ：,将以上方程写成矩阵形式就是：,预测反褶积,处理模块,在Omega系统中，预测反褶积处理由以下三个模块完成，即：,1、预测反褶积谱分析（PRD_DCN_SPCTRL_ANL) ；,2、预测反褶积算子设计(PER_DCN_OPR_DESIGN)：,3、反褶积算子应用(DCN_OPR_APPLY)。,预测反褶积处理模块 在Omega系统中，预测反褶,输入数据,预测反褶积谱分析,预测反褶积算子设计,反褶积算子应用,叠加,输入数据,输入数据预测反褶积谱分析 预测反褶积算子设计 反褶积算子应用,预测反褶积谱分析,功 能：计算给定时窗的自相关函数。该自相关函数用于构造计算预测滤波算,子的方程：,主要参数：1、确定时窗 的参数（起始时间、时窗长度）：,根据资料情况和处理目的确定。,2、自相关长度：,可根据算子长度确定：自相关半长度不能小于算子长度，因,为方程右端用到了自相关函数的值，所以自相关半长度至少,要等于算子长度加预测距离。,预测反褶积谱分析,预测反褶积算子设计,功 能：计算预测反褶积算子。每个道的各个时窗都有自己的预测反褶积,算子。如果需要，你也可以先将谱分析输出的自相关函数按某种,方式（如炮集）进行叠加，然后设计统一的算子。,输 入：自相关函数文件,输 出：1、预测反褶积算子；2、估算子波=预测反褶积算子的反信号。,主要参数：,1、算子长度：这里指的是预测滤波算子的长度。,预测反褶积算子长度=预测滤波算子+预测距离-1,算子长度确定预测滤波方程组的阶数，而构成预测滤波方程组的数据来源于自相关函数。如果自相关函数半长度小于算子长度，则不足部分补零，而如果自相关函数半长度大于算子长度，则后面多余部分不用。,算子长度到底该怎样选择？,在反射系数为白噪声的条件下，记录的自相关函数就是子波的自相关函数，所以最好将整个子波自相关函数用来构建方程组，太短都太长不能正确反映子波的特性。但实际记录的自相关函数并不真等于子波自相关函数，因为反射系数并不是理论上白噪声，记录的自相关函数可能有多个峰值，用以计算预测反褶积算子的自相关函数应尽可能避开第二个峰值。,预测反褶积算子设计,2、预测距离,预测距离即前面提到的是一个重要的参数，它对反褶积的功能起决定性作用。预测反褶积是输出预测误差：,由该式可以看出，时间t+ 处的预测反褶积结果是子波的前个点与t+ 处前面的个反射系数的乘加结果，后面的子波根本未用。这相当于子波被截断成个点再与反射系数函数褶积。,越小，子波被截得越短，反褶积的功能越强，当=1时，预测反褶积变成了脉冲反褶积；反之，反褶积的功能越弱，当大于子波长度时，预测反褶积不起作用。在实际资料处理中，的大小可根据要求达到的频带宽度确定。,但既然预测反褶积是截断子波，对子波就有一的要求，也就是子波能量须集中在前端，这就是为什麽预测反褶积要假定子波是最小相位的缘故。,实际资料的地震子波通常都不理想的最小相位，所以预测距离不可任意给小，否则会改变有效信号的能量和波形特征。,2、预测距离,3、算子修改,本模块设计的算子是最小相位，它不改变输入的相位特征。可修改为,（1）零相位 ：相位谱为0，该算子不改变输入的相位谱,（2）纯相位 ：振幅谱为1，该算子不改变输入的振幅谱,缺省：不改变,3、算子修改,4、白化因子,预测方程并非在任何情况下都可已求解。该方程有唯一确定解的条件是：它的系数矩阵是正定的，即它的各子行列式的值都大于0。由于这里的系数矩阵是自相关函数构成的，所以可以保证它的各子行列式的值都不小于0，即它应该是半正定的。为了使系数矩阵变为正定，以便求解方程，我们就将矩阵的对角线元素增加一个百分数B，将预测方程改造为：,以上做法实际上是将x(t)的自相关函数加一个能量为B的脉冲函数，这相当于在地震记录x(t)上加一个白噪声，故称这一改造为预先白噪化。B称为白噪系数或白噪因子。,4、白化因子,在实际应用中，并非仅仅是为了使方程有唯一确定解。因为数学上的解有时并不适合实际物理问题的要求。以脉冲反褶积为例：,在脉冲反褶积中，反褶积算子a(t)是地震子波b(t)的反信号：,b(t)*a(t)=,(t),在频率域就是,：,或,显然要使上式成立，对任何频率f,必须有B(f) 0，此外， B(f) 也不太接近0，否则会使A(f)的值在这一频率点上过大。因此预先白噪化在频率域就是将B(f) 加上一个小 数w，使其不那么接近0，这时有：,W就是白噪因子。,在实际应用中，并非仅仅是为了使方程有唯一确定解。,在反褶积处理中，大多都需要预先白噪化这一步骤，但它不是反褶积理论推导中的必然步骤，而是根据我们的需要人为地加上去的。因此白噪因子只能是一个很小的数，即只能对理论公式做少许修改。过大的白噪因子可能把理论公式改得面目全非。,在实际处理中，白噪因子参数一般应根据资料的具体情况由试验确定，但通常在参数允许的范围内给一个较小的数即可。,在反褶积处理中，大多都需要预先白噪化这,反褶积算子应用,功 能： 用业已设计好的反褶积算子完成反褶积处理。你用某,个模块计算出了某种类型反褶积算子，都可用该模块完成,应的反褶积运算。例如，当反褶积算子来源于地表一致性,反褶积算子设计时，该模块将完成在地表一致性反褶积处,理。,输入文件：,1、地震数据文件 PRIMARY_SEIMIC,2、反褶积算子文件 SECONDARY_OPERATORS,输出文件：,1、反褶积后的地震输出 FILTERED_SEISMIC,2、未做反褶积的地震道输出 UNFILTERED_SEISMIC,3、算子输出 OPERATORS,反褶积算子应用,地表一致性反褶积,基本原理和计算方法,Omega系统中地表一致性反褶积的实现,地表一致性反褶积基本原理和计算方法,基本原理和计算方法,Omega系统的地表一致性反褶积的反褶积的算法与预测反褶积相同，但获取构建预测方程自相关函数的方法不同，它采用地表一致性算法计算子波的自相关函数，基本原理和计算方法如下：,地震记录x(t)可表示为子波w(t)与地层反射系数函数y(t)的褶积再加上噪声n(t)：,由于地表的不一致性，各道记录的子波w(t)并不一样。记j点激发i点接收的子波为w,ij,(t),则有：,这里：s,j,(t) = 带有炮点影响的子波分量：激发条件对子波的滤波作用,r,i,(t) = 带有检波点影响的子波分量：接收条件对子波的滤波作用,g(t),(i+j)/2,= 与共中心点有关的子波分量：反射点的地质因素对子波的,滤波作用,m(t),(i-j)/2,= 与偏移距有关的子波分量：偏移距（入射角）对子波的滤,波作用,反褶积学习基础课件,对上式做付氏变换，并略去下标后，得到子波复频谱的表达式：,W（f）= S(f) R(f)G(f) M(f),将各个频谱写成指数形式：,W(f)= A,w,e,i,w,S(f)= A,s,e,i,s,R(f)=A,r,e,i,r,G(f)= A,g,e,i,g,M(f)= A,m,e,i,m,其中，字母A表示振幅谱，字母表示相位谱。,这样一来，子波复频谱的振幅谱和相位谱可写为：,A,w,=A,s,A,r,A,g,A,m,(振幅谱),w,=,s,+,r,+,g,+,m,(相位谱),对振幅谱两端取对数，得：,LnA,w,=LnA,s,+LnA,r,+LnA,g,+LnA,m,（对数振幅谱）,对上式做付氏变换，并略去下标后，得到子波复频谱的表达式：,当炮点和检波点变化时，对数振幅谱方程将变成一系列方程。方程个数一般大于未知数个数，但其中独立的方程个数一般小于未知数个数，因此该方程组无确定解。为了求解该方程组，用最小平方法，即要求输入谱与分解谱的误差能量达到最小。于是可得出四个分量的迭代计算公式，再用Gauss-Seidel 法计算出子波对数振幅谱的炮点、检波点、共中心点和共偏移距分量。,当炮点和检波点变化时，对数振幅谱方程将变成一系列方,高斯-赛德尔迭代法,为书写方便，将对数振幅谱表达式中的符号稍加改变。,记：A,w,=A, A,s,=S, A,r,=R, A,g,=G, A,m,=M,将各分量分解出来，要求分解后各分量之和与原来对数振幅谱的误差总能量,达到最小。,由此得到以下四个高斯-赛德尔迭代方程：,其中，i是共炮点序号，j为共检波点序号，,k=(i+j)/2 为共检波点序号，l=i-j为共偏移距序号。,高斯-赛德尔迭代法,将一道记录子波的各分量加在一起，该道子波的对数振幅谱就得到了。作为地表一致性反褶积，应把要去掉的分量加在一块，得到要去除滤波因素的对数振幅谱，再取幂、平方后做反付氏变换得到其自相关函数。有了自相关函数就可以构造预测滤波方程，解此方程 求出预测滤波算子进而得到反褶积算子。,如果我们将四个分量都加起来，那么得到的是整个地震子波w(t)的自相关函数，用它求得的反褶积算子能够压缩子波，起到提高分辨率的作用。当预测距离参数置为一个采样间隔，其结果就是地表一致性意义下的脉冲反褶积。这是地表一致性反褶积的第一层意义，即使用地表一致性反算法的反褶积。,如果我们只将炮点和检波点分量加起来，那么得到的是炮点和检波点分量和的自相关函数，用它求得的反褶积算子能够压缩炮点和检波点的滤波函数，起到减小炮点和检波点的滤波的作用。当预测距离参数置为一个采样间隔，其结果是将炮点和检波点的滤波作用变为单位脉冲，炮点和检波点的滤波作用被完全去掉。这是地表一致性反褶积的第一层意义，即用反褶积方法实现子波的地表一致性。,将一道记录子波的各分量加在一起，该道子波的对数振幅谱就,Omega系统中地表一致性反褶积的实现,在OMEGA系统中，地表一致性反褶积处理由四个模块组成：,1、地表一致性反褶积分析,2、地表一致性反褶积谱分解,3、地表一致性反褶积算子设计,4、反褶积算子应用,地表一致性反褶积分析,地表一致性反褶积谱分解,地表一致性反褶积算子设计,反褶积算子应用,输入地震数据,输入地震数据,地表一致性反褶积处理流程图,Omega系统中地表一致性反褶积的实现 在OMEG,一 地表一致性反褶积分析,功 能：计算输入地震数据指定时窗数据的对数功率谱。本模块采用自相关,或自回归谱分析（最大熵谱分析）方法计算功率谱，对数功率谱的,开平方得到对数振幅谱。,由自相关函数求功率谱：记录,功率谱=自相关函数频谱,由自回归谱分析,求功率谱,：,自回归谱分析也称最大熵谱分析。一个p阶的自回归过程X(n)可表示为：,其中，,a(k) 是自回归过程的模型参数或称回归系数；,U(n) 是一个白噪扰乱序列，它的自相关函数为脉冲函数。,（1）式表明，X(n)的当前值可用其前p个值的线性组合再加上U(n)的当前值来表示。,自回归过程X(n)的功率谱为：,这里，,一 地表一致性反褶积分析,(2)式中的为扰动序列U(n) 的方差，（3）式中的T为过程序列的采样间隔。,在自回归谱分析中，自回归模型参数a(k)可以由过程数据估算出来。于是可由（2）式和（3）式求出过程的功率谱。 过程X(n) 的自相关函数r(n) 可表示为：,上式两端做付氏变换：,自相关的频谱就是过程本身的功率谱：,自相关r(n)也是一个p阶自回归过程，即：,由（4）、（5）两式也可分别计算功率谱和模型参数。,(2)式中的为扰动序列U(n) 的方差，（3）式中的,令A(n)=1,a(n),h(n)为A(n)的逆，则（1）式可写为,如果x(n)是地震数据，则(6)式正好是地震数据的褶积模型描述。h(k)是,地震子波，u(k)是反射系数序列。而A(n)是地震子波的逆。,根据以上自回归谱分析的原理，模块提供了两种计算功率谱的方法：,Yule-Walker 法和inverse of the inverse法。,Yule-Walker 法的实现步骤如下：,1）计算自相关函数,2）计算自回归模型参数：地震子波的逆A(n)即是脉冲反褶积算子，于是用脉冲反褶积求A(n),，,从而得到自回归模型参数a(k),3）根据（2）式和（3）式求功率谱,注：自相关函数的最大延迟值等于自回归阶数p,inverse of the inverse法的实现步骤如下：,1）用与Yule-Walker 法相同的方法计算A(n),2）求子波:子波=A(n)的最小平方逆,3）计算功率谱: 功率谱=子波自相关的付氏变换,令A(n)=1,a(n),h(n)为A(n)的逆,输 入：叠前地震记录。,输 出：各道指定时窗的对数振幅谱,主要参数：,1、地表一致性反褶积参数组(*SURF_CON_DECON),（1）固定自相关半长（ CONSTANT_ACOR_LENTH,）,本参数为所有时窗指定一个自相关计算长度（从0时移到最大时,移），它将自动把已经给定的各时窗的自相关长度参数覆盖。如果,该参数小于时窗长度，将自动改为时窗长度。,大于0的数，或,USE_WNDW_ACOR_LENTH：,用USE_WNDW_ACOR_LENTH参数组指定,缺省：无,（2）自回归模式(AUTOREGRESSIVE_MODEL),这个参数指定谱分析的方式。,YULEWALK 用Yulewalk法,INVERSE 用inverse of the inverse法,NO_MODELING 不用自回归模式（直接用自相关函数）,省略值： Yulewalk,输 入：叠前地震记录。,（3）白噪因子(PRE_WHITE_PCT),这个白噪因子用在自回归谱分析，而非地表一致性反褶积所用的白噪因子，那个白噪因子在地表一致性反褶积算子设计中给出。所以如果不用自回归模式，则本参数可勿略。,缺省：0.01,（4）最小时窗长度(MIN_WNDW_LENGTH),本参数指定计算自相关函数的最小时窗长度值,如果你在时窗设计参数组中给的时窗短于本参数,时窗将自动增加到这个长度。在某些位置如记录尾部,时窗不能加长到这个长度时,该时窗的自相关函数就不计算。这个参数通常不需给出，只有当你的处理要求非常长的自相关长度时才需要。,取值范围： 大于等于固定自相关函数半长度,缺省：最小自相关长度的4倍（COMPUTED,）,（3）白噪因子(PRE_WHITE_PCT),2、时窗自相关长度参数组(*WNDW_ACOR_LENGTHS),当多个时窗的自相关长度各不相同时，则用列表方式给出不同时窗的自相关长度。下面是一个表的例子：,a,2,称a是最小相位延迟信号,2、若a,1,a,2,称a是最大相位延迟信号,3、若a,1,=a,2,称a是等延迟信号,任一n+1项信号 b=(b,0,b,1,b,n,)可分解为n个两项信号 的褶积。,如果,1、所有两项信号 都是最小相位延迟信号，则b是最小相位,2、所有两项信号 都是最大相位延迟信号，则b是最大相位,3、既有最大相位延迟也有最小相位延迟，则b是混合相位,信号的相位特征也可用其z变换来定义：,1、 z 变换的根都在单位圆外，信号是最小相位,2、 z 变换的根都在单位圆内，信号是最大相位,3、单位圆内外都有根，信号是混合相位,最小相位信号的其它定义：,以上有关最小相位信号的定义数学意味比较浓，听起来比较抽象，下面介绍最小相位信号的几个性质，对处理者或许有用。,1、物理可实现的最小相位信号的反信号也是物理可实,现的最小相位信号；,2、最小相位信号与其它信号者褶积，其结果的相位特,征与被褶积的原信号相同；,3、最小相位信号的能量主要集中在时间的前端。,以上有关最小相位信号的定义数学意味比较浓，听起来,调谐反褶积(HARMONIZER_DECON),功 能,Omega系统的调谐反褶积是一种叠后处理模块，它的功能有两个：,（1）提高地震道的纵向分辨率；,（2）衰减多次波。本模块可衰减二次和三次反射类型的多次波，不考虑三次以上反射的多次波。,本模块所使用的计算反褶积算子的方法为预测反褶积。与常规的预测反褶积不同的是：调谐反褶积对各道的每一个样点都分别设计反褶积算子。这就做到了完全的时空变，因而特别适用于地震子波有时-空变化的情况。,显然，用一个样点是无法设计反褶积算子的。事实上，它必然要用包含这个样点在内的一个时窗（通常以该点为时窗中点），只不过时窗向下移动时，每次只移动一个