第2章 反褶积-1.doc

第二章 反褶积 反褶积是借助压缩基本地震子波来改善时间分辨率的一种处理过程。为搞清这一过程要求综合研究正演问题，即必须首先研究记录的地震道的积木式分段单元。地层是由不同类型岩性的岩层组成的，每种岩石类型都有地球物理学家所可利用的某种物理特性。至于地震勘探，则根据波传播速度和岩层密度确定岩层。密度与速度的乘积称之为地震波阻抗，地震资料分析期望的最终成果就是地震波阻抗剖面。我们有在井中直接检测岩层速度和密度的方法，这种方法能向我们提供地震波阻抗与深度的关系。在地面上沿测线记录到的地震反射波就是由于两地层之间的波阻抗差引起的。记录到的反射记录可通过反射率与震源子波的褶积来模拟。 下面分别对褶积模型、各种反滤波进行介绍，并给出应用实例。2.1 褶积模型 我们从图1给出的一个实际声测井记录入手，该声测井曲线是层速度与深度的关系图。实际的速度测量是以 2英尺的采样间隔在10005400英尺之间的深度段内完成的。借助简单的斜坡把速度函数外延至地面。该声测井记录显示出明显突变和强低频趋势特征，这两者构成了总的速度变化。实际上我们通常用CMP道集作速度分析进行估算的就是这种低频趋势。对声测井曲线可通过人工分段提取其速度趋势，其结果可列表如下： 由声测井记录确定的层速度趋势 表1地层序号层速度(ft/s)深度范围(ft)12100010002000219000 20002250318750225025004126502500377551965037755400 实际上该层速度是逐渐减小的。 我们所做的就是形成一组恒定层速度的层组。把测井曲线进行这种分段多少有点类似于地质家对假想的地下模型所做的分层。地质家是根据岩性分层，而我们根据声测井曲线的分段性质提取的分层则是以速度差为依据的。下面对表1中所确定的地层的岩性分类：地层序号岩 性12345灰 岩泥质灰岩(泥岩含量逐渐增加)泥质灰岩泥 岩白 云 岩 在声测井曲线的低频趋势上附加有高频分量。这些急剧变化的波动表明岩层局部特性的变化属性。例如灰岩层可能含有砂泥岩夹层。岩层中的孔隙度变化可能产生较低的层速度。我们还要注意到，测井资料确实精度有限，有些高频变化，特别是以周波跳跃的形式出现的，并不是真实的。不过在实际工作中，我们尽力消除这些问题。 速度和密度测井资料为地球物理特征转换成地下实际地质剖面提供了依据。现在我们的工作是确定测井记录与实际的地震记录道的关系。地震波阻抗定义为密度和速度的乘积。根据实际观测，我们了解到，与速度随深度的变化相比，密度的变化是非常小的。因此，我们可以假定波阻抗差基本上与岩层间速度差是相同的。当然也未必都是这种情况，但这种假设在目前的讨论中，使我们可以只用声测井记录。建立地震道正演模型的第一组假设是：假设l a地层由许多常速水平层组成； b震源产生的平面纵波，沿法线入射到地层界面，在这种情况下不产生横波。 在构造复杂地区，第一种假设是不成立的。第二种假设要求我们必须用零炮检距资料工作，而实际上我们从来都不用自激自收的方法采集数据。但是另一方面，如果地层界面较深，相对来说排列长度较小，我们就可以认为在已知界面上的入射角很小，基本上可以忽略反射系数对角度的依赖。根据假设1a、b，我们定义与界面有关的反射系数，比如说第1层和第2层之间的反射系数r=I(2)I(l)/I(2)+I(l)，式中I是与每层有关的由密度和速度乘积给定的地震波阻抗。这里我们还假设密度随深度是不变的，因此，e=v(2)-v(1)/v(2)+v(1)。在垂直入射时，反射系数是反射波振福与入射波振幅之比。进而，从它的定义我们看出，反射系数就是声阻抗变化与两倍的平均声阻抗之比。如果v(2)大于v(1)，则反射系数为正；如果v(2)小于v(1)则反射系数为负。在图1中还绘出了由声测井曲线v(2)计算出的反射系数e(z)(z是深度变量)。每个尖脉冲位置给出了地层界面的深度，而每个尖脉冲的大小则给出界面反射波振幅占入射的单位振幅平面波的振幅的比例。 现在我们把由声测井记录求得的反射系数序列转换成时间序列。我们选择一个采样间隔比如说2毫秒，那么知道了速度函数之后，其深度轴就可以转换成双程垂直时间轴，如图1(c)所示。反射时间序列e(t)的常规显示也用变面积加波形(相同的序列重复6次以增加强反射的识别度)的形式显示出来。该反射系数序列e(t)代表由等时间间隔、即采样率所划分的一系列假想地层界面的反射率，这种地层模型称之为Goupilland模型。因为地震波是按时间记录的，所以我们需要用这种地层模型工作。由e(t)可以识别主要的反射，由地层2和3之间的界面产生的反射约在0.3秒，而来自4至5层之间界面的反射约在0.5秒。 反射系数序列e(t)仅由“一次”反射组成。我们必须把层外及层间多次反射波考虑在内以获取更真实的反射率序列。如果我们有一个单位振幅的尖脉冲震源，那么记录到的地震记录将代表地层的脉冲响应，该脉冲响应包含着一次反射和多次反射。应用Kunetz方法获得的脉冲响应示于图1(d)，该图同时给出了其变面积加波形显示。 现在我们考虑诸如炸药或气枪这样的脉冲震源的情况。在一个无限小的时间间隔之后，震源时间函数变成了有一定延续时间的子波，并且是有限带宽的。图2给出了一个实例。当这个波向地下传播时会发生两种效应。第一是由于球面发散效应而使整个振幅衰减，其衰减速率与1/r成正比，这里r是震源到地下某一点之间距离。第二是使高频衰减的岩层吸收效应，在前面(第一章)我们曾经研究过这种现象。在时间深度域中基本地震子波的演变如图2所示，所以在任意给定时间上的子波都与震源起始激发时的子波不同。波形随时间的这种变化称之为非平稳性。如何考虑这个问题呢？我们试图用多时窗设计反褶积因子来解决这种非平稳性。另一种方法是作时变的谱白噪化。这些在后面将给予详细讨论，下面作出第二个假设：假设2：震源波形在地下传播时是不变的，即是平稳的。 参见图3。我们从图3(a)显示的震源的脉冲响应开始，该子波沿深度方向传播，在双程时间0.2秒时遇一地层界面，该界面的反射系数用图3(b)中的一个尖脉冲表示。该子波自身重复，并由反射系数标定其大小和极性，如果我们有由图3(d)-图3 (f)的各个尖脉冲所表示的一系列地层界面的话，那么该子波便以相同的方式作自身重复。如果反射系数是负的，该子波会以相反的极性作自身重复，如图3(t)所示。现在来考虑这些反射率的总体效应，如图3(g)所示。这组稀疏尖脉冲序列对基本子波的响应可将各个脉冲响应直接叠加在一起而获得。这个线性过程称之为“叠加原理”，它是由基本子波与图3(g)所示反射率序列的褶积获得的，在第一章我们已经用数值例子证明了褶积过程，在图3(g)中显示的稀疏尖脉冲序列对基本子波的响应具有一些重要特征。值得指出的是，对0.2秒和0.3秒处的尖脉冲，能识别出二个界面的存在，但根据这个综合响应却不能识别出三个相距很近的反射界面。为了区分出这些界面，我们必须恢复震源子波波形来获得稀疏的尖脉冲序列，这正好是我们可据以获得反射率对基本子波的响应的褶积的反过程。因此，非常合理的把这个反过程称之为“反褶积”。 现在对由图1(d)的声测井记录得到的脉冲响应应用叠加原理。代表震源一时间函数的基本子波与脉冲响应相褶积产生合成地震记录，图4对这一过程作了描述，图1(e)也显示了这样做的结果。我们求出的这种一维零炮检距地震记录不含随机环境噪声。实际的地震记录应该加上噪声，如图4所示。现在我们已经完善了对地震记录的“褶积模型”。这种“褶积模型”的积木式结构由下面方程描述： X(t)W(t)*e(t)n(t) (2.1)式中：X(t)：地震记录，W(t)：基本震源子波 e(t)：反射系数序列，n(t)：随机环境噪声方程(2.1)是图4描述的褶积模型的基本数学表达式。我们希望用反褶积方法确定的就是这种反射系数序列(严格地说，是脉冲响应)。 地震记录上的随机噪声有几种来源。外部源包括风和环境噪音，如在检波器附近的牲畜和车辆的通过和行走，或检波器埋置的不好。噪声的内部源来自记录仪器本身。图5是纯噪声地震记录和它的特征。严格地讲，随机信号具有一白谱，它包含着所有的频率成分。因此，这就意味着自相关函数在零延迟时为一个尖脉冲，而在所有其他延迟时为零。由图5可以看到，这些特征要求得到了很好的满足。你可能永远得不到完善的随机信号，因为它若是完善的就势必是以类似于该图的形式出现。 现在我们进一步研究褶积模型方程。在这个方程中已知仅是地震记录X(t)，我们需要估算反射率e(t)。我们不知道震源子波W(t)。在某种情况下我们可以记录震源子波，例如：实际工作中对空气枪组合的子波是要进行记录的。但是，这仅仅是在震源组合激发一开始时的波形。我们已经知道子波形状随着它向地层深部传播而发生变化，所以实际上我们研究的是一种有效子波，它代表总体的、多少是平均的震源-时间函数。最后还有我们事先了解不到的环境噪声n(t)。这样，我们就有了三个未知量、一个已知量的一次方程。我们能解这个方程吗？有人认为不能，但在实际上，对地震资料应用反褶积可看作是常规处理的一个重要步骤，同时多年来我们发现它确实是提高时间分辨率的一种有效方法。那么对未知的e(t)又如何求解呢？让我们进一步作些假设： 假设3：噪声分量n(t)是零。 假设4：震源子波是已知的。 现在我们可以求解反射系数序列e(t)，因为我们有一个方程和一个未知量。在下一节，我们研究求e(t)的一种方法。频率域褶积模型 无噪声地震记录由褶积模型表示成下列方程： x(t)w(t)*e(t) (2.2)时间域的褶积相当于频率域的乘积。准确地说，地震记录的振幅谱等于震源子波的振幅谱与反射系数谱的乘积，图6是脉冲响应e(t)、震源子波w(t)和地震记录x(t)的振幅谱。子波的振幅谱与地震记录的振幅谱之间的相似性是相当明显的。在数学上，这意味着除了差一个比例因子外，x(f)=w(f)，也就是说，假定R(f)为比例因子。在研究了反射系数序列的振幅谱以后，实际上我们可以看到它们分布在整个频谱带宽上，并围绕一个常数值在上、下波动。正如在图5中见到的那样，在整个频谱带宽范围内具有平直谱的时间序列代表一个随机过程。但对于代表地层响应的反射率序列来说，这并不总是正确的。例如旋回沉积在反射率中引入周期性分量。但由声测井记录获得的大量的反射率观测值表明，一般说来。地层反射率的随机性是其突出的普遍特征。 把这个讨论引伸到震源子波、反射率和地震记录的自相关函数。由图6可见，地震记录和基本子波的自相关函数相当类似。当然这种类似性限于子波自相关不为零的那件延迟时间。实际上，反射率的随机性表明，除了零延迟时外，它的自相关对于所有的延迟时都为零。图5显示的随机嗓音也是这种情况，但在随机序列的自相关和反射率的自相关之间存在一个细微的差别。比较图5和图6可以看出，反射率的自相关总是在极大值的零延迟之后跟着有一个相当大的负延迟值，而随机噪音就没有这种情况。除此之外，两个自相关的总的特征相似是相当明显的。根据我们对频谱和自相关函数所作的这些观测结果，可以作出下面假设： 假设5：反射率是一个随机过程，意指在它们的自相关和振幅谱类似的这一点上地震记录具有震源子波的特征。 在这一节后面我们会发现这个假设对实现预测反褶积是一个关键问题。因为它允许我们用已知的地震记录代替未知的震源子波的自相关。因此，有了假设5就不需要假设4了，也就是说，为了从地震记录中提取脉冲响应我们不必知道震源子被。我们感兴趣的是子波的逆，而且我们发现维纳滤波仅要求自相关函数。根据假设5，对于该函数我们可以采用地震记录的自相关。2.2 反滤波 确定一个滤波因子a(t)，使它与已知的地震记录x(t)的褶积产生e(t)的估算值，这个滤波过程可用下式描述； e(t) a(t)*x(t) (2.3)将(2.3)式代人(2.2)式，则 x(t)w(t)*a(t)*x(t) (2.4)方程两边消去x(t)，得到下面的表达式： (t)w(t)*a(t) (2.5)这里(t)为克罗内克函数； (t)= (2.6)我们可以解方程(2.5)求出滤波因子a(t) a(t)=(t)*w(t) (2.7)这里w(t)是已知震源子波w(t)的逆(反子波)。所以，为了由地震记录计算出反射率，我们要求出的滤波因子，在数学上变成了震源子波的反子波。方程2.5的意义在于当反滤波用于基本子波时，将其转换成t=0处的尖脉冲。同样，对地震记录应用反滤波时，将其转换成确定反射率的一系列尖脉冲。所以，如果已知震源子波，那么，反滤波就是一种反褶积方法。下面将反滤波概括于表1。 如何求震源子波的反子波呢？在数学上，是应用Z-变换来完成的。现在我们举一个简单的数值例子。假定基本子波是两点时间序列：(1，-1/2)。那么该子波的Z-变换就可由下面的多项式确定； w(Z)=1-(1/2)Z (2.8)变量Z的幂是与时间序列中每个采样点有关的单位时间延迟的个数。第一项为零延迟，则“Z”的幂为零，第二项为 l个单位的延迟，则“Z”的幂为 1。所以 Z变换是“Z”的多项式，其系数正是原来时间序列的元素。对于初次见到Z变换时的人来说，这似乎是神秘的。实际上，Z-变换与傅里叶变换之间有着自然的关系。Z变量被定义成Z=exp(idt)，其中是角频率，dt是采样间隔，子波的逆w(t)由对(2.8)式给出的Z变换进行多项式除法获得。 W(Z)=l/(1-Z)=1+Z+Z2+ (2.9)这个反时间序列是W(Z)的系数，即(1，1/2，1/4，)，这就是滤波因子a(t)值得注意的是，虽然它们衰减的很快，但它具有无限多个系数。和任何一种滤波处理一样，要截断滤波因子，首先我们考虑二点滤波因子(l，l/2)，把它与子波褶积，如表2所示： 反滤波因子(1，12)与震源子波(1，12)褶积 表2 0 1 -1/2 输出 1/2 1 0 1 0 l/2 1 0 0 1/2 1 1/4理想的结果是(1，0，0)，而实际输出是(1，0，-1/4)。虽然这不是最期望的结果，但比输入子波(1，1/2)确实尖锐了。 让我们看看通过对反滤波增加一项系数能否使效果更好些。由三项滤波得到的输出如表3所示： 反滤波因子(1，1/2，1/4)与震源子波(1，l/2)的褶积 表3 1 -1/2 输出 -1/4 1/2 1 1 -1/4 l/2 1 0 -1/4 l/2 1 0 -1/4 1/2 1 -1/8三点滤波的输出是(1，0，0，-1/8)，与二点滤波输出(1，0，-1/4)相比，(1，0，0，-1/8)是理想输出(1，0，0)的更精确的表示。至此，我们可以说，在三点滤波输出的情况下，有较少的能量漏失到非零延迟上，所以输出是更加脉冲化的。当把更多的项包括到反滤波因子中时，其输出就变得更加接近于在延迟处的脉冲，然而实际上，因子的长度总有一个点数的限制，这样其输出永远不会是一个完善的尖峰脉冲。 在输入子波为(1，-1/2)的情况下，其反子波具有迅速衰减至零的系数。那么当输入子波为(-1/2，1)时会怎样呢？这时它的反滤波因子是由一个发散序列(-2，-4，-8)给定的。要注意，该系数随着时间向后延迟迅速增大。现在我们截断这个序列并用二点因子与输入子波(-1/2，1)作褶积： 反滤波(一2，4)与输入子波(12， 1)的褶积 表4 -1/2 1 输出 -4 -2 1 -4 -2 0 -4 -2 -4实际输出是(1，0，4)，而期望输出是(1，0，0)。显然这一结果不仅与期望输出相差甚远，而且还不如输入子波(-1/2，1)尖锐。造成输出结果不好的原因是：反滤波系数随时间延迟增大而不是衰减；而在截断时，我们实际上是截掉了较大的系数。2.3 最小平方反滤波 当有一个特征良好的输入子波时，用上节描述的反滤波能产生较好的接近于尖脉冲的输出。现在的问题是，我们能否做得更好一点。让我们来研究下面的问题：假设输入子波为(l，1/2)，我们要找到一个二项滤波因子(a，b)，使其实际输出与期望输出(1，0，0)之间的误差是在最小平方意义下的最小值。将滤波器(a，b)与输入子波(1，-1/2)褶积计算实际输出。 滤波器(a，b)与输入子波(1，-1/2)的褶积 表5 1 -1/2 实际输出 期望输出 b a a 1 b a b-1/2a 0 b a -b/2 0误差的累加能量定义为实际输出与期望输出系数之差的平方和： E=(a-1)2+(b-a/2)2+(-b/2)2 (2.10)我们的任务是寻求系数(a，b)，使E取最小值，这就要求E关于系数(a，b)的变分趋于零，简化并取E对“a”和“b”的偏导数并使其等于零，得到： (5/2)a-b=2 (2.11a) (5/2)b-a=0 (2.11b) 这两个方程有两个未知数，即滤波系数(a，b)，可以把这个方程组写成简便的矩阵形式。 (2.12)对此方程求解滤波系数，我们得到：(a，b)=(20/21，8/21)将该滤波器用于输入子波，其结果如下表： 输入子波(1，-1/2)与滤波系数(20/21，8/21)的褶积 表6 1 -1/2 输出 8/21 20/21 20/21 8/21 20/21 -2/21 8/21 20/21 -4/21 为了对该结果的脉冲化程度作定量评价，并与表2列出的反滤波结果进行比较，我们计算两者的误差能量。 二项滤波和最小平方滤波的误差 表7 输 入：(1，-1/2) 期望输出：(1，0，0) 实际输出 误差能量 反滤波 -1 0 -1/4 1/16 最小平方滤波 20/21 -2/21 -4/21 1/21 在力图将输入子波(1，-1/2)整形成在零延迟处的尖脉冲(1，0，0)的过程中，最小平方滤波能产生较小的误差。 我们再研究一下在输入子波为(-1/2，1)时最小平方滤波的性能。还记得，对于这个子波(表4)反滤波产生了不稳定的结果。我们要找到二项滤波因子(a，b)，使它与输入子波(-1/2，1)褶积后产生期望脉冲输出(1，0，0)的估计值，并希望实际输出与期望输出之间的平方误差为最小。 滤波因子(a，b)与输入子波(-1/2，l)的褶积 表8 -1/2 1 实际输出 期望输出 b a -a/2 1 b a -b/2+a 0 b a b 0其误差的累积能量为： E=(-a/2-1)2+(-b/2+a)2+b 2.13整理且取E对“a”和“b”的偏导数并令其等于零，我们得到： 5/2a-b=-1 (2.14a) 5/2b-a=0 (2.14b)写成矩阵形式： 2.15解出滤波系数，(a，b)(-10/21，-4/21)，将该滤波器用于输入子波得出如下结果： 输入子波(-1/2，1)与滤波系数(-10/21，-4/21)褶积 表9 -1/2 1 输出 -4/21 -10/21 5/21 -4/21 -10/21 -8/21 -4/21 -10/21 -4/21 和第一个子波一样，我们也对反滤波和最小平方滤波的结果作一定量分析。 二项反滤波和最小平方滤波的误差 表10 输 入：(-1/2，1) 期望输出：(1，0，0) 实际能量 误差能量 反滤波 -1 0 -4 16 最小平方滤波 5/21 -8/12 -4/21 336/441最小平方滤波还是比反滤波误差小得多，但是当输入子波为(-1/2，1)时，两种滤波器都产生较大误差。下一节讨论其原因。2.4 最小相位 在2.2节和2.3节我们曾用了两个输入子波分别对反滤波和最小平方反滤波进行了数值分析。子波1是1，(1/2)，子波2是(1/2)，1，分析结果表明，子波1转换成零相位脉冲所产生的误差比用子波2时的误差小，子波1比子波2更接近于零延迟脉冲(l，0，0)，也就是说，子波2比子波1更接近于延迟脉冲(0，l，0)。所以，从直观上讲，我们可以得出这样的结论：如果我们选择的期望输出非常类似于输入序列的能量分布的话，那么其误差将减少。子波1的能量主要集中在波始；相反，子波2的能量则集中在波尾。由第一章可以知道，通过直接改变相位延迟谱来改变子波形状而不用修改振幅谱。 如果子波的能量主要集中在波始，我们称该子波是最小相位；同样，如果子波的能量集中在波尾，我们说它是最大相位；如果子波的能量主要在中心部位，我们称该子波是混合相位。我们发现把子波定义成具有一定寿命的瞬时波形是有意义的。最小相位子波是单边的，在t=0之前它等于零，这与我们的直观感觉是一致的，即系统对某个激发的响应是产生在激发之后，我们把这种t0子波为零的子波称为“因果关系”。所以最小相位子波是符合因果律的，因此它具有有限寿命，并是可实现的。 让我们对这些观测作定量分析，现在考虑下面4个3点子波： 子波A：(4， 0，-1) 子波B：(2，3，-2) 子波C：(-2，3，2) 子波D：(-1，0，4)计算某一时刻每个子波的累加能量： 子波 不同时间采样的能量 0 1 2 A 16 16 17 B 4 13 17 C 4 13 17 D 1 1 17 图2.4.1是这些数值曲线图，所有4个子波具有相同的总能量，即17个单位。这说明它们的振幅谱是相同的，但是能量增长的速率对于每个子波是不同的。例如，子波A，在最初的时间延迟其能量便迅速地逼近它的总能量值。对于子波B和C的能量增长比较缓慢。对于子波D的情况，其能量以最慢的速率聚集。由图2.4.1我们见到，对于子波A和D的能量曲线形成了上、下边界。我们发现子波A具有最小能量延迟。相反，子波D具有最大的能量延迟，这就是罗滨逊(Robinson)能量延迟理论。具有最小能量延迟的子波称为最小延迟子波；具有最大能量延迟的子波称为最大延迟子波。 我们知道，时间延迟相当于相位延迟，所以我们研究这些子波的相位延迟谱是相当重要的。图2.4.2是相位延迟谱。由图可见，子波A通过频率轴具有最小的频率变化，因此，我们说子波A是最小相位。子波D具有最大的相位变化，因此它是最大相位。最后是子波B和D，它们的相位变化界于两者之间，我们说它们是混合相位。 因为所有这四个子波具有相同的振幅，所以它们的功率谱亦相同，它们也应该具有相同的自相关，通过下表的数据可以证明这一点。可见自相关的零延迟等于每个子波包含的总能量，即17个单位对于任何子波，功率谱下的面积都等于自相关函数的零延迟值，帕斯韦尔(Parseval)定义严格地阐述了这一点。 子波A、B、C、D的自相关输出 表1 4 0 -1 输出 4 0 -1 17 4 0 -1 0 4 0 -1 -4 2 3 -2 2 3 -2 17 2 3 -2 0 2 4 -2 -4 -2 3 2 -2 3 2 17 -2 3 2 0 -2 3 2 -4 -1 0 4 -1 0 4 17 -1 0 4 0 -1 0 4 -4 正如在2.2和2.3节所见到的那样，脉冲反褶积是把震源子液压缩成零延迟脉冲的过程。其精度不仅取决于滤波器的长度(与2-和3-点滤波相比较)，而且还取决于子波是否是最小相位。严格地讲，脉冲反褶积因子是子波的逆。如果子波是最小相位，我们能得到一个稳定的反子波，它也是最小相位的。所谓稳定，是指滤波器系数形成一个收敛的序列。如子波(1，-1/2)的情况，它的逆为(1，1/4，1/2，)。这个逆表示一个稳定的脉冲反褶积滤波器。另一方面，如果子波是最大相位，我们就不能得到稳定的逆。子波(-1/2，l)就是这种情况，由发散序列(-2，-4，-8，)给出它的逆。最后，如果子波为混合相位，我们仍然得不到稳定的逆。这个讨论引出下面一种假设：假设6：震源子波是最小相位的，所以它有最小相位的逆。脉冲反褶积要求作这种假设。下面我们将看到，子波即使不是最小相位的，它也可能转换成延迟脉冲。2.5 最佳维纳滤波 让我们取第一类的期望输出(1，O，0)，这在研究反滤波和最小平方滤波 时已经考虑过了，我们计算输入子波(1，12)的自相关： 输入子波(1，-12)的自相关 表2 1 -1/2 输出 1 -1/2 5/4 1 -1/2 -1/2 该输出与方程(2.12)左边的22矩阵中的第1列相同，比例因子为2。换句话说，方程(2.12)可改写成如下形式： (2.16)两边除以2可得到： (2.17) 现在来计算期望输出(1，0，0)与输入子波(1，12)的互相关。 期望输出(1，0，0)与输入子波(1，1/2)的互相关 表3 1 0 0 输出 1 -1/2 1 1 -1/2 0 该输出与方程(2.17)右边的这列矩阵相同。总之，左边的矩阵元素是输入子波互相关延迟,而右边的矩阵元素是期望输出与输入子波的互相关延迟。 让我们对子波完成类似的运算。计算它的自相关： 输入子波(-1/2，1)的自相关 表4 -1/2 1 输出 -1/2 1 5/4 -1/2 1 -1/2 我们已经注意到子波(-1/2，1)的自相关与子波(1，-1/2)的自相关是相等的(表2.3.7)。正如在第4节已经讨论的那样，这是具有如同振幅谱的一组子波的重要特征。矩阵方程2.15重写成： (2.18)两边除以2得到： (2.19) 现在我们来计算期望输出(1，0，0)与输入子波(-1/2，1)的互相关。 期望输出(1，0，0)与输入子波(1，1/2)的互相关 表5 1 0 0 输出 -1/2 0 -1/2 -1/2 1 0可以看到其输出与方程(2.19)右边的相同。方程左边的矩阵元素是输入子波的自相关延迟。为导出最小平方反滤波我们用丁矩阵方程(2.17)和(2.19)，然后把这些滤波器用于输入子波，其目的是把输入子波压缩成零延迟脉冲。我们发现，方程(2.17)和方程(2.19)左边的矩阵是由输入子波的自相关延迟组成。另外，右边的这列矩阵是由期望输出(即零相位脉冲与输入子波)的互相关延迟组成。实际上，为了得到把输入转换成“任意”期望的输出的滤波器，诺贝尔特维纳推广了这些观测值。滤波器长度“n”的矩阵方程的一般形式为： (2.21)式中r(i)，a(i)和g(i)，i=0，1，2，3(n1)分别为输入子波的自相关，滤波器系数和期望输出与输入子波的互相关。这种滤波器a(i)称之为最佳维纳滤波器，而最佳是指实际输出与期望输出之间的平方误差为最小的意义上来讲的。我们已经注意到，当期望输出是零延迟脉冲(1，0，00)时，维纳滤波相当于最小平方“反滤波”，换言之，后者是前者的一个特例。维纳滤波器不仅限于零延迟脉冲，而且可用于研究任意期望输出的大量问题。我们通过研究期望输出是零延迟脉冲的特例，已经得到了标准方程的一般形式(2.20)。左面的流程图概括地描述了维纳滤波器的设计和应用。期望输出可能的实际选择： 类型1：零延迟脉冲， 类型2：任意延迟的脉冲， 类型3：输入序列的时间超前形式， 类型4：零相位子波， 类型5：某种期望的任意形状， 对维纳滤波器系数的标准方程组求解需要将左边的矩阵进行反转。我们能通过手算解22矩阵的数学例子。然而对于长滤波器因子来说，矩阵求逆变得非常昂贵。由方程(2.2O)可见，自相关矩阵是对称的，这种特殊矩阵称之为托布里兹矩阵。我们使用莱文森递归图表来求解滤波器系数。实际上我们先计算出二点滤波器，由此导出三点滤波器，以此类推，直到导出n点滤波器为止。莱文森递归是非常快的，尤其是应用了阵列处理机以后。实践中，以最佳维纳滤波器理论为基础的滤波算法就是著名的维纳-莱文森理论。脉冲反褶积 获得类型1的期望输出的处理过程称为脉冲反褶积。期望输出(1，0，O，O)与输入子波x(0)，x(1)，x(2)，x(n1)的互相关产生了x(0，0，0，O)。取标准方程组(220)的一般形式的特殊式： (2.21) 方程(221)由1/x(0)定标。在第3节讨论的最小平方反滤波具有与上面矩阵方程相同的形式，所以，脉冲反褶积等同于最小平方反滤波。 图2.5.1是用维纳-莱文森算法对脉冲反褶积的概括。我们从混合相位子波(a)幅开始，其振幅谱(b)幅表明，该子波是带限的，大部分能量限制在1050Hz的范围内。相位谱是跳跃的，表明它不是最小相位谱(C幅)。我们用方程2.21的自相关函数计算脉冲反褶积因子(e幅)。该因子的振幅谱就是子波振幅谱的逆(f幅)。我们期望，脉冲反褶积的目的是通过加强弱频率分量使谱拉平。该因子的相位谱证明它是最小相位(g幅)，把这个因子用于输入子波，我们得到(k)幅。理想情况下，我们能得到(a)幅那样的零相位脉冲。问题出现在哪里呢？其答案在于假设6。我们是由混合相位输入子波开始的。在(h)幅上所显示的是反褶积因子的逆。实际上，这是(a)幅上显示的混合相位子波的最小相位等效值。它们两个都具有相同的振幅谱(b幅和i幅)但它们的相位谱确有很大差别(C幅和j幅)。因为脉冲反褶积相当于反滤波，所以最小相位子波是反褶积因子的逆。因此，该因子的振幅谱是最小相位子波振幅谱的逆(f幅和i幅)，并且前者的相位谱是后者的负值，如g幅和j幅所示。如果震源子波是最小相位，那么提取该震源子波的方法就是直接计算脉冲反褶积因子并取其逆。总之，如果我们没有最小相位子波，我们就不能把它转换成纯零延迟脉冲(k幅)，尽管振幅谱实际上是平的(l幅)，但输出的相位谱并不是最小相位的(m幅)。最后，应该记住，脉冲反褶积因子是非最小相位输入子波的最小相位等效值。预白噪化 我们已经知道，脉冲反褶积因子的振幅谱是输入子波振幅谱的逆。图2.5.2对此作了描述。在振幅谱上如果有些零点那将如何呢？实际上，我们取最小相位子波并应用带通滤波(图2.5.3)，可见到，滤波后的子波及褶积不产生纯尖脉冲，而是得到为高频的前趋波和后续波所伴随的脉冲。显然，这个结果是不能用的，因为反褶积因子试图增加不存在的频率，如同输入振幅谱所见到的那样。这种问题在地震记录上可能产生吗？我们很少在输入振幅谱上有“空白”的情况。在地震记录上总有噪音存在，在时间域和频率域它都是外加的。此外，产生的数值噪声在频率域也是附加的，但是，为保证数值的稳定性，在反指积之前，我们可以人为地选择“白噪”能级，这称之为白噪化，图2.5.2对此作了描述。在输入子波的谱上加上一定量级的白噪相当于绝大多数加到了自相关的零延迟上。如果白噪化的百分比由数“e”给定，那么所求解的标准方程组(据方程2.21)可能为下面形式： (2.22) 后面我们将研究白噪化能级(e)的选择对反褶积特性的影响。维纳整形滤波 我们已经了解到把子波(-1/2，1)压缩成零延迟脉冲(l，0，0)，其脉冲反褶积会有困难。根据能量分布，该子波更类似于(0，l，0)的延迟脉冲，而不同于零延迟脉冲(1，0，0)，所以，设想把子波(-1/2，1)转换成延迟脉冲的滤波器应产生比零延迟脉冲情况更小的误差是合理的。按照表2.3.11的流程图给出的步骤设计一个滤波器并用来将(-1/2，1)转换成延迟脉冲(0，1，0)。 期里输出(0，1，0)与输入子波(1/2，1)的互相关 表6 0 1 0 输出 -1/2 1 1 -1/2 1 -1/2 由表2.3.9我们已经知道了输入子波的自相关。将表2.3.9和2.3.12的结果代人矩阵方程(2.2O)得到： (2.23)求解该滤波器系数，我们得到(a，b)(16/21，-2/21)，将该滤波器用于输入子波： 输入