反褶积_学习基础.ppt

欢迎各位光临 地震资料处理中的反褶积处理 黄大云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 之间的距离越大 分辨率越高 反之分辨率越低 通常震源产生的信号 震源子波 是较短的 但它在传播过程中 因为大地的滤波作用会逐渐拉长 以至分辨率越来越低 为了提高分辨率 只有两种办法 加大g t 之间的距离或者压缩b t 的延续长度 g t 之间的距离是客观存在 显然我们无法也不应该去改变它 为了提高地震记录的分辨率 只有压缩地震子波b t 的长度 理想的情况是将b t 缩为单位脉冲函数 这时地震记录x t 就是反射系数序列g t 把b t 缩为单位脉冲函数的方法通常是用某种办法设计出一个算子 它与b t 褶积的结果就是单位脉冲函数 由此我们得反褶积的定义 4 反褶积的一般定义反褶积就是去掉地震记录中大地的滤波作用的一种处理方法 所以反褶积也叫反滤波 它用的运算方法归根到底仍然是褶积 但现在的反褶积已不局限于去除大地的滤波作用 凡是对地震子波进行改造的处理都叫它反褶积 5 反褶积处理的目的提高地震记录的分辨率是反褶积处理的目的之一 但对叠前反褶积而言 它却不是主要目的 叠前反褶积的主要目的是使地震子波波形一致 以便获得好的叠加效果 未做反褶积 地表一致性反褶积后 预测反褶积的原理和计算方法 主要内容什么叫预测预测的条件预测滤波预测反褶积预测反褶积的计算Omega系统的预测反褶积模块及其主要参数 因为Omega系统的几个重要的反褶积模块 如预测反褶积 地表一致性反褶积和调谐反褶积都用的是预测算法 所以有必要先说一说预测反褶积的原理和计算方法 什么叫预测 预测就是根据过去和现在已发生的事实判定将来会出现的情况 在数学上 对一个时间函数的预测是指该函数某一点的值用其前面若干个值的线性组合表示出来 这种预测称为线性预测 预测的条件并非所有事物都可线性预测 函数x t 可线性预测的条件是 x t 为平稳随机过程 即它的统计特征 数学期望和方差是与时间无关的量 且自相关函数rxx 只与时差 有关 我们认为地震记录满足以上条件 因而可做预测 预测滤波 在地震勘探中 我们认为地震记录是平稳随机过程 因而可以预测 根据地震记录褶积模型的假设 地震记录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 为一物理可实现的最小相位信号 因此有 当t 0时 a t 0将g t a t x t 带入x t 得 令s j k 上式变为 再令得到 上式为一褶积表达式 它说明 x t 是c s 对x t 的过去和现在值的滤波结果 称它为x t 的预测值 c s 称为预测滤波因子 实际值与预测值的差e t x t x t 称为预测误差 叫做预测间隙 预测步长或预测距离 预测反褶积 因为 所以预测误差为 特别地 当 1时 有 e t 1 b 0 g t 1 该式表明 当预测距离等于1时 预测误差与反射系数只差一个常数因子 因而可视为反射系数 于是 只要在预测滤波中输出预测误差就达到预测反褶积的目的 这就是预测反褶积 但通常不用 1这种理想情形 而是令 为大于1的某个数 当 1时 预测反褶积就是脉冲反褶积 预测反褶积的计算 预测反褶积计算的关键是求得预测滤波因子c s 由于子波未知 不能用公式直接计算 可用最小平方法 最小平方法的数学模如下 输入信号 x t 设预测滤波因子 c s c 0 c 1 c m 期望输出 x t 0 预测输出 预测误差 误差总能量 选取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 预测反褶积谱分析 功能 计算给定时窗的自相关函数 该自相关函数用于构造计算预测滤波算子的方程 主要参数 1 确定时窗的参数 起始时间 时窗长度 根据资料情况和处理目的确定 2 自相关长度 可根据算子长度确定 自相关半长度不能小于算子长度 因为方程右端用到了自相关函数的值 所以自相关半长度至少要等于算子长度加预测距离 预测反褶积算子设计功能 计算预测反褶积算子 每个道的各个时窗都有自己的预测反褶积算子 如果需要 你也可以先将谱分析输出的自相关函数按某种方式 如炮集 进行叠加 然后设计统一的算子 输入 自相关函数文件输出 1 预测反褶积算子 2 估算子波 预测反褶积算子的反信号 主要参数 1 算子长度 这里指的是预测滤波算子的长度 预测反褶积算子长度 预测滤波算子 预测距离 1算子长度确定预测滤波方程组的阶数 而构成预测滤波方程组的数据来源于自相关函数 如果自相关函数半长度小于算子长度 则不足部分补零 而如果自相关函数半长度大于算子长度 则后面多余部分不用 算子长度到底该怎样选择 在反射系数为白噪声的条件下 记录的自相关函数就是子波的自相关函数 所以最好将整个子波自相关函数用来构建方程组 太短都太长不能正确反映子波的特性 但实际记录的自相关函数并不真等于子波自相关函数 因为反射系数并不是理论上白噪声 记录的自相关函数可能有多个峰值 用以计算预测反褶积算子的自相关函数应尽可能避开第二个峰值 2 预测距离预测距离即前面提到的 是一个重要的参数 它对反褶积的功能起决定性作用 预测反褶积是输出预测误差 由该式可以看出 时间t 处的预测反褶积结果是子波的前 个点与t 处前面的 个反射系数的乘加结果 后面的子波根本未用 这相当于子波被截断成 个点再与反射系数函数褶积 越小 子波被截得越短 反褶积的功能越强 当 1时 预测反褶积变成了脉冲反褶积 反之 反褶积的功能越弱 当 大于子波长度时 预测反褶积不起作用 在实际资料处理中 的大小可根据要求达到的频带宽度确定 但既然预测反褶积是截断子波 对子波就有一的要求 也就是子波能量须集中在前端 这就是为什麽预测反褶积要假定子波是最小相位的缘故 实际资料的地震子波通常都不理想的最小相位 所以预测距离不可任意给小 否则会改变有效信号的能量和波形特征 3 算子修改本模块设计的算子是最小相位 它不改变输入的相位特征 可修改为 1 零相位 相位谱为0 该算子不改变输入的相位谱 2 纯相位 振幅谱为1 该算子不改变输入的振幅谱缺省 不改变 4 白化因子预测方程并非在任何情况下都可已求解 该方程有唯一确定解的条件是 它的系数矩阵是正定的 即它的各子行列式的值都大于0 由于这里的系数矩阵是自相关函数构成的 所以可以保证它的各子行列式的值都不小于0 即它应该是半正定的 为了使系数矩阵变为正定 以便求解方程 我们就将矩阵的对角线元素增加一个百分数B 将预测方程改造为 以上做法实际上是将x t 的自相关函数加一个能量为B的脉冲函数 这相当于在地震记录x t 上加一个白噪声 故称这一改造为预先白噪化 B称为白噪系数或白噪因子 在实际应用中 并非仅仅是为了使方程有唯一确定解 因为数学上的解有时并不适合实际物理问题的要求 以脉冲反褶积为例 在脉冲反褶积中 反褶积算子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 SEIMIC2 反褶积算子文件SECONDARY OPERATORS输出文件 1 反褶积后的地震输出FILTERED SEISMIC2 未做反褶积的地震道输出UNFILTERED SEISMIC3 算子输出OPERATORS 地表一致性反褶积 基本原理和计算方法Omega系统中地表一致性反褶积的实现 基本原理和计算方法Omega系统的地表一致性反褶积的反褶积的算法与预测反褶积相同 但获取构建预测方程自相关函数的方法不同 它采用地表一致性算法计算子波的自相关函数 基本原理和计算方法如下 地震记录x t 可表示为子波w t 与地层反射系数函数y t 的褶积再加上噪声n t 由于地表的不一致性 各道记录的子波w t 并不一样 记j点激发i点接收的子波为wij t 则有 这里 sj t 带有炮点影响的子波分量 激发条件对子波的滤波作用ri t 带有检波点影响的子波分量 接收条件对子波的滤波作用g t i j 2 与共中心点有关的子波分量 反射点的地质因素对子波的滤波作用m t i j 2 与偏移距有关的子波分量 偏移距 入射角 对子波的滤波作用 对上式做付氏变换 并略去下标后 得到子波复频谱的表达式 W f S f R f G f M f 将各个频谱写成指数形式 W f Awei wS f Asei sR f Arei rG f Agei gM f Amei m其中 字母A表示振幅谱 字母 表示相位谱 这样一来 子波复频谱的振幅谱和相位谱可写为 Aw AsArAgAm 振幅谱 w s r g m 相位谱 对振幅谱两端取对数 得 LnAw LnAs LnAr LnAg LnAm 对数振幅谱 当炮点和检波点变化时 对数振幅谱方程将变成一系列方程 方程个数一般大于未知数个数 但其中独立的方程个数一般小于未知数个数 因此该方程组无确定解 为了求解该方程组 用最小平方法 即要求输入谱与分解谱的误差能量达到最小 于是可得出四个分量的迭代计算公式 再用Gauss Seidel法计算出子波对数振幅谱的炮点 检波点 共中心点和共偏移距分量 高斯 赛德尔迭代法为书写方便 将对数振幅谱表达式中的符号稍加改变 记 Aw A As S Ar R Ag G Am M将各分量分解出来 要求分解后各分量之和与原来对数振幅谱的误差总能量达到最小 由此得到以下四个高斯 赛德尔迭代方程 其中 i是共炮点序号 j为共检波点序号 k i j 2为共检波点序号 l i j为共偏移距序号 将一道记录子波的各分量加在一起 该道子波的对数振幅谱就得到了 作为地表一致性反褶积 应把要去掉的分量加在一块 得到要去除滤波因素的对数振幅谱 再取幂 平方后做反付氏变换得到其自相关函数 有了自相关函数就可以构造预测滤波方程 解此方程求出预测滤波算子进而得到反褶积算子 如果我们将四个分量都加起来 那么得到的是整个地震子波w t 的自相关函数 用它求得的反褶积算子能够压缩子波 起到提高分辨率的作用 当预测距离参数置为一个采样间隔 其结果就是地表一致性意义下的脉冲反褶积 这是地表一致性反褶积的第一层意义 即使用地表一致性反算法的反褶积 如果我们只将炮点和检波点分量加起来 那么得到的是炮点和检波点分量和的自相关函数 用它求得的反褶积算子能够压缩炮点和检波点的滤波函数 起到减小炮点和检波点的滤波的作用 当预测距离参数置为一个采样间隔 其结果是将炮点和检波点的滤波作用变为单位脉冲 炮点和检波点的滤波作用被完全去掉 这是地表一致性反褶积的第一层意义 即用反褶积方法实现子波的地表一致性 Omega系统中地表一致性反褶积的实现 在OMEGA系统中 地表一致性反褶积处理由四个模块组成 1 地表一致性反褶积分析2 地表一致性反褶积谱分解3 地表一致性反褶积算子设计4 反褶积算子应用 地表一致性反褶积分析 地表一致性反褶积谱分解 地表一致性反褶积算子设计 反褶积算子应用 输入地震数据 输入地震数据 地表一致性反褶积处理流程图 一地表一致性反褶积分析功能 计算输入地震数据指定时窗数据的对数功率谱 本模块采用自相关或自回归谱分析 最大熵谱分析 方法计算功率谱 对数功率谱的开平方得到对数振幅谱 由自相关函数求功率谱 记录功率谱 自相关函数频谱由自回归谱分析求功率谱 自回归谱分析也称最大熵谱分析 一个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 两式也可分别计算功率谱和模型参数 令A n 1 a n h n 为A n 的逆 则 1 式可写为如果x n 是地震数据 则 6 式正好是地震数据的褶积模型描述 h k 是地震子波 u k 是反射系数序列 而A n 是地震子波的逆 根据以上自回归谱分析的原理 模块提供了两种计算功率谱的方法 Yule Walker法和inverseoftheinverse法 Yule Walker法的实现步骤如下 1 计算自相关函数2 计算自回归模型参数 地震子波的逆A n 即是脉冲反褶积算子 于是用脉冲反褶积求A n 从而得到自回归模型参数a k 3 根据 2 式和 3 式求功率谱注 自相关函数的最大延迟值等于自回归阶数pinverseoftheinverse法的实现步骤如下 1 用与Yule Walker法相同的方法计算A n 2 求子波 子波 A n 的最小平方逆3 计算功率谱 功率谱 子波自相关的付氏变换 输入 叠前地震记录 输出 各道指定时窗的对数振幅谱主要参数 1 地表一致性反褶积参数组 SURF CON DECON 1 固定自相关半长 CONSTANT ACOR LENTH 本参数为所有时窗指定一个自相关计算长度 从0时移到最大时移 它将自动把已经给定的各时窗的自相关长度参数覆盖 如果该参数小于时窗长度 将自动改为时窗长度 大于0的数 或USE WNDW ACOR LENTH 用USE WNDW ACOR LENTH参数组指定缺省 无 2 自回归模式 AUTOREGRESSIVE MODEL 这个参数指定谱分析的方式 YULEWALK用Yulewalk法INVERSE用inverseoftheinverse法NO MODELING不用自回归模式 直接用自相关函数 省略值 Yulewalk 3 白噪因子 PRE WHITE PCT 这个白噪因子用在自回归谱分析 而非地表一致性反褶积所用的白噪因子 那个白噪因子在地表一致性反褶积算子设计中给出 所以如果不用自回归模式 则本参数可勿略 缺省 0 01 4 最小时窗长度 MIN WNDW LENGTH 本参数指定计算自相关函数的最小时窗长度值 如果你在时窗设计参数组中给的时窗短于本参数 时窗将自动增加到这个长度 在某些位置如记录尾部 时窗不能加长到这个长度时 该时窗的自相关函数就不计算 这个参数通常不需给出 只有当你的处理要求非常长的自相关长度时才需要 取值范围 大于等于固定自相关函数半长度缺省 最小自相关长度的4倍 COMPUTED 2 时窗自相关长度参数组 WNDW ACOR LENGTHS 当多个时窗的自相关长度各不相同时 则用列表方式给出不同时窗的自相关长度 下面是一个表的例子 WNDW NUMACOR LENGTH1120220031404120 1 时窗序号 WNDW NUM 2 自相关半长度 ACOR LENGTH 自相关一半的长度 0到最大的正延迟 二 地表一致性反褶积谱分解功能 将地表一致性谱分析模块所计算出来的对数功率谱分解为共炮点 共检波点 共中心点和共偏移距四个分量 计算方法 高斯 赛德尔迭代法 输入 前步处理输出的对数功率谱 输出 1 对数功率谱的分量文件 2 分量合并文件主要参数 1 一般参数组 GENERAL 1 迭代次数 MAX NUM ITER 用Gauss Seidel方法求解的迭代次数 缺略值 6 2 第几次迭代前要做共中心点分量平滑 LAST ITER CMP SMOOTH 在有的情况下 方程的解不能很好的收敛 从而出现畸变 对此可用对共中心点分量做横向平滑来克服 本参数确定哪次的迭代结果要做平滑 它一旦指定 从第一次迭代开始到迭代次数等于这个参数时都要做平滑 小于迭代次数的整数NO SMOOTH 不做平滑缺省值 COMPUTED 迭代次数减1 3 CMP的最小重复度 MIN CMP MULT 该参数用来计算横向平滑滤波器的长度 如果该参数值为m 记录的最大道号 MAX TRACE NUM 为n 则平滑滤波器的长度为k n 2m 滤波器的数值如下 1 k 1 k 1 k当k为奇数时 1 2k 1 k 1 k 1 2k当k为偶数时 取值范围 大于等于3缺省值 3 4 检波点间距 GROUP INTERVAL 本参数确定偏移距的分组 如这个参数为N 则偏移距为0到N 2的分为第一组 偏移距为N 2到3N 2的分为第二组 如此等等 当使用缺省值时 它将使用道头中TRC SPACING 道距 值或CMP SPACING 共中心点距 值乘2 如果道头中这两个值均不存在 而又选定了要做共偏移距分量分解 则必须指定该参数 缺省值 COMPUTED 2 分量次序参数组 COMPORDER 本参数组指定参加分解的分量及其次序 在2维处理中 四个分量都应指定 在3维处理中 因不能做共中心点分量平滑 所以不使用共中心点分量 1 炮点分量在Gauss Seidel解中的次序 ORDER SOURCE COMP 该参数确定要否计算炮点分量 如果要 那么它是第几个解 EXLUDE不计算或1 2 3 4中的一个数缺省值 1 2 检波点分量在Gauss Seidel解中的次序 ORDER DETECTOR COMP 该参数确定要否计算检波点分量 如果要 那么它是第几个解 EXLUDE不计算或1 2 3 4中的一个数缺省值 2 3 中心点分量在Gauss Seidel解中的次序 ORDER MIDPOINT COMP 该参数确定要否计算中心点分量 如果要 那么它是第几个解 EXLUDE不计算或1 2 3 4中的一个数缺省值 3 4 偏移距分量在Gauss Seidel解中的次序 ORDER OFFSET COMP 该参数确定要否计算偏移距分量 如果要 那么它是第几个解 EXLUDE不计算或1 2 3 4中的一个数缺省值 4 3 结果输出参数组 COMP OUTPUT 你可以要求输出各分量之和 以对一地震道计算一个反褶积算子 为此 你可在输出栏目中选择SUMMED SPECTRA即可 若在输出栏目中选择个分量的谱 DECOMP SPECTRA 则将根据你的要求输出各个分量 已计算的分量可以不输出 没有计算的分量不能要求输出 1 输出炮点分量 OUTPUT SOURCE COMP YES 输出炮点分量NO不输出炮点分量缺省值 YES 2 输出检波点分量 OUTPUT DETECTOR COMP YES输出检波点分量NO不输出检波点分量缺省值 YES 3 输出中心点分量 OUTPUT MIDPOINT COMP YES输出中心点分量NO不输出中心点分量缺省值 YES 4 输出偏移距点分量 OUTPUT OFFET COMP YES输出偏移距分量NO不输出偏移距分量缺省值 YES 三 地表一致性反褶积算子设计功能 生成地表一致性反褶积算子 方法 用前一模块输出的分量和求出自相关函数 再构造预测滤波方程 解此方程得到预测滤波算子 并根据预测距离参数构成反褶积算子 设预测滤波算子为 a1 a2 an 预测距离为t 则反褶积算子为 1 0 0 0 a1 a2 an 其中0的个数为t 1 输入 对数功率谱分量或分量和 输出 地表一致性反褶积算子文件 估算子波 EstimatedWavelets 这里输出的子波实际上是反褶积算子的逆 即它与反褶积算子的褶积为单位脉冲 所以只有预测距离为一个采样间隔时 它才是子波真正的估算 主要参数 1 预测距离2 算子长度3 算子修改类型它们的意义与预测反褶积相同 4 地表一致性能量校正振幅水平 POWER CORRECTION AMP 这是地表一致性反褶积算子设计参数组 SURF CON OPER 中的一个参数 它确定是否要使反褶积输出道的振幅水平保持等于给定的常数N 如果要 则算子将被扩大或缩小 以使反褶积后数据的均方根振幅水平等于N 能量校正只适用于单个时窗反褶积 如果地表一致性谱分析是多个时窗 而你又要求地表一致性能量校正 则作业将被终止 NO SCPC不做地表一致性能量校正其它给定能量输出水平N缺省 NO SCPC 4 反褶积算子应用 功能 用业已设计好的反褶积算子对输入地震记录做褶积 完成反褶积处理 该模块并非地表一致性反褶积专用 预测反褶积也用的是这一模块 输入 1 地震记录2 反褶积算子输出 反褶积后的地震记录 预测反褶积与地表一致性反褶积的比较相同点 都采用预测算法求反褶积算子不同点 预测反褶积直接用本道的自相关函数作为子波的自相关函数 可在某个域 如炮域 内叠加自相关函数以求该炮域的一个统一的反褶积算子地表一致性反褶积采用地表一致性分解法求各道子波的自相关函数 子波整形反褶积处理模块 基本原理和计算方法Omega系统中子波整形反褶积的实现 基本原理和计算方法子波整形反褶积 以下简称子波反褶积 对震源子波整形 也可用于实现地表一致性的反褶积 它不考虑激发和接收条件对子波具体有什么改造 反正把各个震源子波都变成一个已知波形就行 使用本方法有一个前提 即假定同一炮记录中 各接收点对子波的影响是随机的 同一接收点道集中 各炮点对子波的影响也是随机的 子波反褶积既可在共炮点域进行也可在共接收点域进行 如果分别在两个域进行 则地表一致性的效果会更好 子波反褶积的计算可分为两步 1 子波估算 2 子波整形 一子波估算以共炮点域为例 根据地震记录的褶积模型 地震记录可表示为地震子波与反射系数函数的褶积 在反射系数函数为白噪序列的假设下 子波的功率谱可由记录的功率谱估算 当我们只考虑记录前端的资料的情况下 可将地震子波看成就是震源子波 同一炮记录各道的震源子波是一样的 将各道的功率谱平均即得到统计意义下的震源子波功率谱 功率谱的平方根即是振幅谱 再假定震源子波的相位特征为最小相位 这个最小相位谱可用其对数振幅谱的希尔伯特变换求得 由此得到一个估算的震源子波 估算震源子波的过程如下图 道记录 震源子波 希尔伯特变换 二子波整形震源子波整形有两种方式 将震源子波整形为尖脉冲 把已求的最小相位谱反号 振幅谱取倒数后构成复频谱 再做反付氏变换即得到震源子波整形反褶积算子 2 将震源子波整形为一个期望的目标子波 这时 整形算子的频谱为OP w 可表示为 其中 分子为期望目标子波的复频谱 分母为估算的震源子波复频谱 最后用求得的整形算子对各道记录做褶积 共接收点域的子波反褶积将接收条件变成一致 它的做法与共炮点域相同 因为子波反褶积假定震源子波为最小相位 所以如果要在两个域做子波反褶积 那么第一域的期望目标子波必须指定最小相位 Omega系统中子波整形反褶积的实现 子波反褶积由四个模块共同完成 它们依次是 1 子波反褶积谱分析 DESIG SPCTRL ANL 2 叠加 STACK 3 子波反褶积算子设计 DESIG OPR DESIGN 4 反褶积算子应用 DCN OPR APPLY 子波反褶积谱分析对地震炮记录作谱估算 叠加模块按共炮集 或共接收点集 对谱估算进行叠加求其均值 子波反褶积算子设计用谱估算的平均值设计一个反褶积算子 最后反褶积算子应用模块用这个反褶积算子与该炮的各道做褶积 一子波反褶积谱分析 DESIG SPCTRL ANL 功能 本模块对地震炮记录作谱估算 计算方法 有三种谱的类型可供选择选择 1 自相关 2 功率谱 3 对数功率谱 以10为底 有两种谱估算方法可供选择 1 直接谱分析 用时窗内的数据直接做付氏变换得出复频谱 并由此计算功率谱或对数功率谱 用直接谱分析不输出自相关 2 间接谱分析 间接谱分析先求时窗数据的自相关 自相关的付氏变换即是功率谱 取对数即得对数功率谱 用间接谱分析可选择以上三种输出类型的任一种 直接谱分析的分辨率较高但统计的稳定性较差 间接谱分具备统计的稳定性但分辨率较低 子波反褶积主要参数本模块的主要参数均在子波反褶积分析参数组 DESIGNATURE DECON 中 共有以下12个参数 1 谱估算方法 SPECTRAL METHOD 本参数指定做谱估算的方法 ACOR INDIRECT基于自相关的间接谱分析DIRECT直接谱分析缺省 ACOR INDIRECT2 输出谱类型 OUTPUT TYPE 本参数指定谱估算输出类型 ACOR自相关POWER功率谱LOG POWER对数功率谱缺省 ACOR3 计算起始时间延迟 START TIME DELAY 本参数用来计算分析时窗的起始时间 道的起始时间加上这里给出的值即为分析时窗的起始时间 选择本参数时 注意避开直达波和折射直达波 缺省 COMPUTED 200毫秒 4 子波反褶积时窗长度 WNDW LENGTH 这个参数决定分析时窗长度 缺省值2000毫秒效果不错 建议使用缺省值即可 缺省 COMPUTED 2000毫秒 5 固定自相关半长度 CONSTANT ACOR LENGTH 这里给出的长度是从0延迟到最大延迟的自相关长度 它等于自相关长度一半加1 如果你给出的值大于时窗长度 模块会自动将其置为时窗长度 取值范围 大于0 小于时窗长度的数缺省 无6 自相关函数镶边类型 ACOR TAPER 对自相关函数做镶边处理可提高谱估算的分辨率和统计的稳定性 这个参数选择镶边函数的类型 NO TAPER不做镶边处理BARTLETTBartlett三角镶边GAUSSIANGaussian函数镶边缺省 NO TAPER 7 最小快速付氏变换长度 这个参数指定付氏变换的长度 样点个数 付氏变换的长度越大 频率域的采样率越小 因而频域分辨率越高 但这会增加计算工作量 因此建议用户在本模块中不使用该参数 而在子波反褶积算子设计模块中使用 当你不指定付氏变换的长度时 模块将根据时窗内的样点数自动确定恰当的变换的长度 本参数仅当你的输出谱为功率谱或对数功率谱时才需要 取值范围 大于固定自相关半长度 采样间隔缺省 NOT USED 不用 8 时窗补0长度 ZERO AMP PAD LENGTH 这个参数指定要在时窗后面加0的个数 它与最小快速付氏变换长度参数有相同的作用 你在时窗后面加了若干0值 增加了时窗长度 自然要增加付氏变换长度 在多数情况下都不必使用这一参数 同样 本参数仅当你的输出谱为功率谱或对数功率谱时才需要 取值范围 大于0的数缺省 NOT USED 不用 9 频率加密采样 FREQ SUPERSAMPLING 为了加密频率采样须增加付氏变换长度 本参数给定一个幂指数 以使原付氏变换长度按本参数的乘方加长 例如 原付氏变换长度为256 本参数的值为3 则新的付氏变换长度为本参数也是仅当你的输出谱为功率谱或对数功率谱时才需要 10 频率平滑带宽 FREQ SMOOTHING 本参数指定对谱做平滑的宽度 平滑公式如下 其中 L 本参数值 单位 赫兹 INPUT 平滑前的谱OUTPUT 平滑后的谱本参数仅用于直接谱分析 缺省 NOT USED 不用 11 时窗前端镶边长度 BEGIN TAPER LENGTH 本参数指定时窗前端镶边长度 镶边类型为余弦函数 取值范围 大于等于0 小于时窗长度的一半缺省 NOT USED 不用 12 时窗后端镶边长度 END TAPER LENGTH 本参数指定时窗后端镶边长度 镶边类型为余弦函数 取值范围 大于等于0 小于时窗长度的一半缺省 NOT USED 不用 前面已说过 镶边可能改变谱分析结果 通常不用为好 二叠加 STACK 如果做共炮点域的子波整形 就对子波反褶积谱分析模块输出的谱分析文件按共炮点叠加 如果做共检波点域的子波整形 就对子波反褶积谱分析模块输出的谱分析文件按共检波叠加 显然 平均值的对数不一定等于对数的平均值 所以对数功率谱不合适叠加 因此在对子波反褶积谱分析中应选择自相关或功率谱类型 三子波反褶积算子设计 功能 子波反褶积算子设计模块用子波反褶积谱分析输出的谱估算经叠加模块做了炮集 或共检波点集 平均的结果为每炮记录 或共检波点道集 设计一个子波反褶积算子输入 1 已叠加的谱估算文件 2 外部目标子波 可选 输出 1 滤波算子 2 道集平均估算子波 主要参数1 子波反褶积算子设计参数组 DESIG OPER 本参数组由以下参数组成 1 白噪因子 PRE WHITE PCT 这个参数表示所估算的子波振幅谱总能量的一个百分比 在计算反褶积算子的过程中 子波振幅谱将用作分母 为避免分母为0 在分母上加一个小的数是必要的 但不能太大 取值范围 大于等于0 小于等于5 0缺省 0 01 2 期望输出相位 OUTPUT PHASE 本参数指定你希望整型后输出子波的相位 在选择这个参数时应考慮后续处理对子波相位特征的要求 MIN最小相位ZERO零相位SPIKE尖脉冲 3 期望输出极性 OUTPUT POLARITY 按照美国地球物理学会的标准 负极性对应于记录检波器的压缩 正极性对应于记录检波器的舒张 这意味着 对于蒸汽枪震源 正反射的初至在野外带上记录为一个负数 西方地球物理公司的转换极性是 正极性表示 正反射系数对应于正数 负反射系数对应于负数 负极性与此相反 这里的极性选择按照西方地球物理公司的转换极性 POSITIVE正极性NEGATIVE负极性缺省 POSITIVE西方地球物理公司的转换极性与美国地球物理学会的标准正好相反 为了使处理结果的极性与野外记录一样 这里宜选用负极性 4 算子预测分量的长度 ANTICIPATION LENGTH 这个参数用来确定算子的零点位置 用算子总长度的百分比表示 算子的零点把算子分为两部份 零点之前称为预测分量 零点之后 包含零点 称为记忆分量 取值范围 大于等于0 小于等于100 缺省 25 2 子波反褶积目标子波参数组 DESIG TARGET 这一组参数也是用来确定目标子波 它可代替仪器参数 如果外部目标子波文件已给出 则此组参数不用 1 目标子波低截频率 TARGET LO CUT 取值范围 大于等于0缺省 仪器响应低截频 2 目标子波低截衰减 TARGET LO SLOPE 以每倍频程分贝方式给出 取值范围 大于等于0 小于180缺省 仪器响应低截衰减 3 目标子波高截频率 TARGET HI CUT 取值范围 大于等于0缺省 仪器响应高截频 4 目标子波高截衰减 TARGET HI SLOPE 以每倍频程分贝方式给出 取值范围 大于等于0 小于180缺省 仪器响应高截衰减3 反褶积算子修改参数组 OPER MODIFY 1 算子修改类型 MODIFY TYPE NO MODIFICATION不修改ZERO PHASE修改为零相位PHASE ONLY修改为纯相位 4 算子归一化开关 OPER NORMALIZE 这个参数指定你是否要对所设计的算子做归一化处理 所谓归一化即算子的每个样点都除以算子的均方根 均方根有两种求法 a 算子的全部样点的均方根 b 只用给定频带范围求均方根 YES要做归一化NO不做归一化缺省 YES5 算子归一化的低截频 NORMALIZE LO FREQ 指定为做归一化求均方根的低端频率 它与算子归一化的高截频参数必须同时使用 取值范围 大于等于0 小于等于算子归一化高截频参数缺省 NOT USED 不用 6 算子归一化的高低截频 NORMALIZE HI FREQ 指定为做归一化求均方根的高端频率 它与算子归一化的低截频参数必须同时使用 取值范围 大于等于0 大于等于算子归一化低截频参数缺省 NOT USED 不用 子波反褶积和地表一致性反褶积比较子波反褶积把所有子波变成相同的目标子波 因而更具一致性 但可能失去原记录子波的特征 地表一致性反褶积 如果目的是子波的地表一致性处理的话 去掉的是地表的不一致对子波的影响 因而原记录反射子波的特征得以保留 但子波的一致性可能不及子波反褶积 前面所述的三种叠前反褶积方法 预测反褶积 地表一致性反褶积和子波反褶积都假设子波为最小相位子波 预测反褶积的公式是在子波为最小相位的前提下推导出来 而在子波反褶积中 用希尔伯特变换估算的是子波的最小相位谱 那么什么是最小相位信号呢 最小相位信号最小相位信号的全称是最小相位延迟信号 可见它与相位延迟有关 1 相位延迟设信号的频谱为 经时间延迟以后变成 则的频谱为 这说明是经滤波器滤波后的结果 该滤波器的相位谱为 由此可得延迟时间与滤波器相位谱的关系 就是由滤波器相谱引起的时间延迟 称此为相位延迟 这是滤波过程为简单时移的情况 对于一般的滤波器 滤波器的相位延迟定义为 该定义说明 滤波引起的时移与频率有关 2 群延迟滤波器的群延迟定义为 群延迟是相位谱函数的微分乘以系数 这说明群延迟描述了相邻频率成分相位谱的变化率 3 最小相位信号的定义一个信号称为最小相位信号 如果其它任何具有相同振幅谱但相位谱不同的信号 它们的群延迟总是大于原信号的群延迟的话 也即是说 在振幅谱确定的情况下 改变相位谱可得若干个 甚至无穷个 信号 其中群延迟最小的那个相位谱就是该振幅谱对应的最小相位谱 用该相位谱得到的信号就是最小相位信号 这样的定义说明最小相位信号的什么性质呢 设最小相位信号为x1 t 它的群延迟为其它具有相同振幅谱的信号为x t 它们的群延迟为于是有或这说明最小相位信号相位谱在任何频率值处的变化率都是最大的 也就是说 构成该信号的各谐波成份的初始相位变化梯度最大 最小相位信号的其它定义 设一两项信号a a1 a2 则1 若a1 a2 称a是最小相位延迟信号2 若a1 a2 称a是最大相位延迟信号3 若a1 a2 称a是等延迟信号任一n 1项信号b b0 b1 bn 可分解为n个两项信号的褶积 如果1 所有两项信号都是最小相位延迟信号 则b是最小相位2 所有两项信号都是最大相位延迟信号 则b是最大相位3 既有最大相位延迟也有最小相位延迟 则b是混合相位信号的相位特征也可用其z变换来定义 1 z变换的根都在单位圆外 信号是最小相位2 z变换的根都在单位圆内 信号是最大相位3 单位圆内外都有根 信号是混合相位 以上有关最小相位信号的定义数学意味比较浓 听起来比较抽象 下面介绍最小相位信号的几个性质 对处理者或许有用 1 物理可实现的最小相位信号的反信号也是物理可实现的最小相位信号 2 最小相位信号与其它信号者褶积 其结果的相位特征与被褶积的原信号相同 3 最小相位信号的能量主要集中在时间的前端 调谐反褶积 HARMONIZER DECON 功能Omega系统的调谐反褶积是一种叠后处理模块 它的功能有两个 1 提高地震道的纵向分辨率 2 衰减多次波 本模块可衰减二次和三次反射类型的多次波 不考虑三次以上反射的多次波 本模块所使用的计算反褶积算子的方法为预测反褶积 与常规的预测反褶积不同的是 调谐反褶积对各道的每一个样点都分别设计反褶积算子 这就做到了完全的时空变 因而特别适用于地震子波有时 空变化的情况 显然 用一个样点是无法设计反褶积算子的 事实上 它必然要用包含这个样点在内的一个时窗 通常以该点为时窗中点 只不过时窗向下移动时 每次只移动一个样点 输出文件 1 SEISMIC一次波反褶积输出 2 AUTO CORR A与第一延迟时相对应的各点的自相关函数 3 AUTO CORR B与第二延迟时相对应的各点的自相关函数 主要参数 1 反褶积算子相位类型反褶积算子一经确定 其相位特征就已确定 该参数无意改变反褶积算子的相位特征 而是指出 你要求按什么样的相位特征来应用反褶积算子 如果认为反褶积算子是最小相位 那么它与地震记录的褶积过程是 1 以前端点为轴将算子褶叠 2 与数据做乘加运算 如果认为反褶积算子是零相位 那就意味着算子是对称的 就按对称算子与数据褶积 2 预测步长时间1求取衰减二次反射多次波的滤波算子的预测步长 从零延迟到延迟为本参数的自相关函数将用来计算滤波算子 它的值一般应取一次反射到第二次反射时间之差 3 预测步长时间2求取衰减第三次反射多次波的滤波算子的预测步长 从零延迟到延迟为本参数的自相关函数将用来计算滤波算子 它的值一般应取一次反射到第三次反射时间之差 4 预测步长时间1的增量 每次迭代时预测步长时间1的增量 5 预测步长时间2的增量 每次迭代时预测步长时间2的增量 6 迭代次数 调谐反褶积可进行迭代运算 即将前次反褶积结果作为输入再做反褶积 谱白化处理 谱白化的基本原理Omega系统中的时变谱白化模块应用谱白化的条件谱白化应在什麽时候使用谱白化后能否保持相对振幅 谱白化的基本原理根据地震记录的褶积模型 地震记录x t 可表示为反射系数函数g t 和子波b t 的褶积 x t g t b t 在频率域就是 X f G f B f 振幅谱关系为 Ax f Ag f Ab f 在记录的振幅谱曲线上 子波的振幅谱表现为低频趋势 反射系数的振幅谱表现为高频变化 子波振幅谱 反射系数振幅谱 记录振幅谱 如果去掉低频趋势 就相当于消除了子波效应 留下的高频变化就是反射系数的振幅谱 这就达到了提高分辨率的目的 通常认为反射系数的振幅谱近乎白谱 所以通过对地震记录振幅谱的白化处理 可得到反射系数的振幅谱 谱白化处理将地震记录振幅谱曲线分成若干小段 设法使每段包围的面积相等 Ax f f 分段越窄 谱越平 但高频细节可能改变 分段越宽 高频细节保持越好 但谱白化效果越差 应用谱白化的条件 谱白化是在假定地震记录的振幅谱等于为反射系数振幅谱与子波振幅谱的乘积 Ax f Ag f Ab f 实际上记录上总有噪声存在 即Ax f Ag f Ab f An f An f 是噪声的振幅谱 因此 谱白化会放大信噪比较低频段的噪声 主要是高低频段的随机噪声 但对能量很强的规则干扰可能有所压制 建议在谱白化前做随机噪声衰减 且加回百分比要小 Omega系统中的时变谱白化模块 功能 将地震记录的振幅谱拉平 使其成为白谱 并可保持原始记录能量随时间的变化特征 实现方法 1 用若干带通滤波器对输入道滤波 2 对每个滤波结果分时窗求振幅包络并平滑以形成增益曲线 3 每个滤波结果与各自的增益曲线相除 4 所有处理后的结果相加 得到该输入道的输出 主要参数 1 带通滤波器个数及其宽度 可选择自动生成和用户给出滤波器 每个滤波器的通频宽度 振幅系数为1的频带范围 不应小于5HZ 每个滤波器的通频宽度越小 输出结果的谱越 白 但过小的振幅宽度 可能抹杀反射系数的振幅谱特征 滤波器的最大振幅宽度越大 谱白化的作用越小 当其大到与输入记录的频带范围相当时 输出与输入就没啥差别了 自动生成滤波器时 最大通频宽度按下式确定 最大通频宽度 高通频 低通频 2 滤波器个数 1 例1 自动生成带通滤波器 FREQ AMP SETAUTO FILT 5 FLIT NUMFREQAMP15 0011011001105 001 例2 用户指定四个带通滤波器 FREQ AMP SETAUTO FILT NO FLIT NUMFREQAMP13 0018113118 001213 00118123128 001323 00128133138 001433 00135140143 001 2 白噪百分比在我们对个滤波道做增益放大时 增益因子是各滤波的振幅包络与输入道振幅包络的比 当输入道振幅包络有0或接近0的值时 会使计算无法进行或者结果不稳定 因此需预先给输入道振幅包络加上一个小的数 就像反褶积一样 即所谓白噪 谱白化应在什麽时候使用 谱白化是一种单道处理 原则上叠前 叠后都可用 但叠前的谱不宜展得太宽 因为子波太短可能影响叠加成像 在应用中发现 叠前使用谱白化展宽了频谱 叠加后频谱又变窄了 这说明 叠前的谱白化提高了噪声水平 经叠加后噪声又被压制了 所以叠前的谱白化有时不能达到预期效果 叠后根据处理要求做谱白化效果要好些 因为叠加道的信噪比一般会大大高于叠前单道的信噪比 谱白化后能否保持相对振幅 根据谱白化的实现方法 它要对每个滤波结果分时窗求振幅包络并平滑以形成增益曲线 每个滤波结果与各自的增益曲线相除得到各频段分量 这样一来 相同时窗内各分量的能量相同 同一分量各时窗的能量也一致 相同时窗内各分量的能量相同 实现了谱白化 同一分量各时窗的能量也一致 就相当于做了时变增益 这必然改变振幅的相对关系 如此看来 谱白化后不能保持相对振幅 但Omega系统的谱白化采用对输出道