第2章 反褶积附录Word下载.docx

第2章 反褶积附录Word下载.docx

《第2章 反褶积附录Word下载.docx》由会员分享，可在线阅读，更多相关《第2章 反褶积附录Word下载.docx（16页珍藏版）》请在冰豆网上搜索。

在已推荐的所有方法中，似乎最小相位反褶积最为突出，其应用最广泛。

尽管最小相位假设本身是关键问题。

我想通过下面的常规处理顺序简单地谈一下当前反褶积问题的进展（Chambers和Larner本人的意见）。

（1）应用几何扩散补偿函数V2t，这将消除由于球面扩散引起的振幅损失。

（2）应用最小相位反Q滤波，这将对衰减进行补偿。

（3）选件。

对海上资料应选用信号反褶积，对可控震源资料应选择将Klauder子波转换成等价的最小相位子波的滤波器。

（4）应用地面一致性反褶积，这将首先着重解决由于在震源和接收点附近的介质非均匀性引起的子波横向变化。

然后解决子波的垂向变化。

（5）叠后和TVF（时变滤波）前应用TVSW。

这将得到一个确切的谱，而没有相位和展平谱对在叠加期间曾被压制的高频部分的影响。

（6）选件。

应用预测反褶积来消除任何短周期的交混回响。

基本准则尽可能多地作确定性的反褶积。

反Q滤波，信号反褶积和把Klauder子波转换为等价于它的最小相位子波的滤波器都是确定性的算子。

我们试图通过统计的方法来解决剩下的问题。

附录：

反褶积的数学基础

A·

1合成地震记录

与地层有关的地震声阻抗定义为I=ρV。

这里的ρ是地层密度，V是地层速度，地震声阻抗的瞬时值由下式给出：

It＝ρtVt（A.1）

考虑平面波理论和垂直入射的情况，界面的反射系数由下式给出：

e=（It＋1-It）/（It+1＋It）（A.2）

如果我们假设密度随深度的变化相对于速度随深度的变化是可以忽略的话，方程（A.2）变为下面形式：

e≌（Vt＋1-Vt）/（Vt+1＋Vt）（A.3）

如果我们有一个单位入射波振幅，那么e（t）的数值与分界面的反射振幅的一部分相对应。

反射系数的时间位置与界面深度成正比。

根据反射系数的知识，我们应用Kunetz方法计算层状地层模型的脉冲响应。

这个脉冲响应不仅仅包括一次反射而且也包括所有可能的多次反射。

最终，脉冲响应与基本震源子波的褶积产生合成地震记录。

我们可以对地震记录再加上随机噪音，但是，在确定反褶积滤波器中使用的褶积模型是没有随机噪音的。

进而，我们假设没有岩层固有衰减。

因此，当震源子波在地层内传播时其波形没有变化。

震源子波同脉冲响应的褶积由下列表达式给出：

X（t）＝w（t）*e（t）（A.4）

X（ω）＝W（ω）E（ω）（A.5）

式中：

X（ω）、E（ω）和W（ω）分别表示地震记录反射率和震源波形的复数傅里叶变换，依据振幅谱和相位谱，这些傅里叶变换表达式如下：

X（ω）＝Ax（ω）eiφx（ω）（A.6a）

W（ω）＝Aw（ω）eiφw（ω）（A.6b）

E（ω）＝Ae（ω）eiφe（ω）（A.6c）

将（A.6）式代入到（A.5）式中，我们有：

Ax（ω）＝Aw（ω）Ae（ω）（A.7a）

φx（ω）＝φw（ω）φe（ω）（A.7b）

这样，作为震源子波与脉冲响应褶积的结果是相位谱相加振幅谱相乘。

我们假设，脉冲响应是“白化处理”，也就是说它的振幅谱是平的：

Ae（ω）＝A0＝常数（A.8）

因此，

Ax（ω）＝A0Aw（ω）（A.9）

这就意味着地震记录的振幅谱是震源子波振幅谱的简单定比形式。

现在我们研究一下x，e和w的自相关函数。

自相关是在不同的时间位置上的时间序列同相轴之间相似程度的一个量度，它是一个求和运算。

其表达式是：

rx（τ）=

，τ=±

0，±

1，…，±

（N-1）（A.10）

这里τ是时间延迟。

随机的时间序列相互之间没有可以相比的同相轴。

rx（τ）＝0τ≠0（A.11a）

rx（0）＝r0＝常数（A.11b）

方程（A.11）说明纯随机序列的自相关除了零延迟之外所有延迟是零。

实际上，零延迟数值是时间序列中累加能量：

r0＝x20+x21+x22+x20+…+x2N-1（A.12）

考虑褶积方程（A.4）的Z变换：

X（z）＝W（z）E（z）（A.13）

我们将1/z放在z的位置并取复数共轭，则有：

（1/z）＝

（1/z）

（1/z）（A.14）

（A.13）和（A.14）等式两边相乘得：

X（z）

（1/z）＝W（z）E（z）

（1/z）（A.15）

等式右边重新整理，则：

（1/z）＝[W（z）

（1/z）][E（z）

（1/z）]（A.16）

最后，定义：

rx=rwre（A.17）

式中rx、re和rw分别是地震记录脉冲响应和震源子波的自相关。

根据“白化”反射率序列的假设，我们有：

rx=r0rw（A.18）

方程（A.18）说明了地震道自相关是震源子波的定比例形式。

而后我们将看到将震源子波转换成零延迟的脉冲需要知道它的自相关。

方程（A.18）说明，即使我们不知道震源子波，也可用地震记录的自相关来代替震源子波的自相关。

A.2震源反子波

反褶积的基本目的是压缩震源波形使其成为一个零延迟时的脉冲，以便于可以分辨密集的反射层，考虑滤波因子f（t），这样

w（t）*f（t）＝σ（t）（A.19）

式中σ（t）是克罗内克（Kronecker）σ函数。

具体地讲，我们可以用震源子波W（t）的倒数来表示f（t）：

f（t）＝1/w（t）（A.20）

我们看到滤波因子f（t）把震源子波转换成一个脉冲，因此，它是一个脉冲反褶积因子。

由于它是震源子波的逆，所以称为“反滤波”。

这个滤波将对地震记录的相似性起作用：

x（t）*f（t）＝w（t）*e（t）*f（t）=e（t）（A.21）

方程（A.21）说明对地震记录使用反滤波时它产生脉冲响应。

我们写出震源子波的Z变换：

W（z）＝W0十W1z十W2z2＋…十Wmzm（A.22）

用多项式相除获得反滤波的z变换：

F（z）=1/W（z）（A.23）

它产生另一个多项式，它们的系数是反滤波的项：

F（z）=f0＋f1z十f2z2＋…十fnzn（A.24）

这里n＋1是滤波长度，多项式除法是一个费时的运算而且滤波长度必须有限。

在震源子波脉冲化的过程中滤波因子的截断产生某些误差。

让我们研究频率域震源反子波。

方程（A.6c）是子波的傅里叶变换，方程（A.19）的傅里叶变换是：

W（ω）F（ω）＝1（A.25）

方程（A.6b）代入上式，我们得到：

F（ω）=1/[Aw（ω）ei

]（A.26）

定义反滤波的傅里叶变换为

F（ω）=Af（ω）ei

]（A.27）

和方程（A.26）比较，有

Af（ω）=1/Aw（ω）（A.28a）

φf（ω）=-φw（ω）（A.28a）

方程（A.28）说明了反滤波的振幅谱是震源子波振幅谱的逆，而它的相位谱是震源子波相位谱的负值。

A.2反滤波

除了方程（A.23）描述的多项式除法之外。

现在我们考虑推导反滤波的另一种不同方法。

震源子波自相关的z变换为：

Rx（z）＝W（z）

（1/z）（A.29）

那么方程（A.19）的z变换就是；

F（z）W（z）＝1（A.30）

由它我们得到：

W（z）＝l/F（z）（A.31）

代人方程（A.29）中：

Rx（z）F（z）＝

（l/z）（A.32）

由于子波是实数时间函数，

rx（t）=rx（-t）（A.33a）

（t）＝w（t）（A.33b）

我们取三点反滤波的特殊情况，方程（A.32）便呈下列形式：

（r2z-2＋r1z-1＋r0＋r1z＋r2z2）（f0＋f1z＋f2z2=W0＋W1z-1＋W2z-2（A.34）

写成简单的矩阵形式，我们有：

（A.35）

这方程左边的方阵中的元素表示未知震源子波的自相关延迟，然而，根据方程（A.18），我们可以用已知的地震记录的自相关延迟来代替。

另外要注意到W（0）是子波在t＝0时的振幅，因此，这里有三个方程四个未知量。

相对于f（0）的归一化，则

（A.36）

这里有三个未知量，即f

（1），f

（2）和V；

和三个方程。

方程（A.36）中的自相关矩阵是一个特殊的矩阵。

首先，它是一个对称矩阵，其次它的主对角线元素是相同的，这叫做托布里兹矩阵。

对滤波系数求解方程（A.36）要求自相关矩阵的逆，对于n个正规方程，其解需要与n2成比例的存贮空间。

而CPU（中央处理机）时间与n3成比例。

可是，由于托布里兹矩阵的特殊性质，莱文森确定的递归方法其要求存贮空间和CPU时间分别与“n”和“n2”成比例。

A.4频率域反褶积

我们从频率域震源子波的自相关开始讨论：

Rw（ω）＝W（ω）

（ω）（A.37）

定义新的函数U（ω）：

U（ω）＝In[Rw（ω）]（A.38）

方程（A.38）两边取幂，则

RW（ω）＝exp[U（ω）]（A.39）

现在再定义另一个函数φ（ω），并重新写方程（A.39）如下：

RW（ω）＝exp{[1/2U（ω）-iφ（ω）][1/2U（ω）+iφ（ω）]}（A.40）

与方程（A.37）相比较可见：

W（ω）＝exp[1/2U（ω）+iφ（ω）]（A.41）

由方程（A.37）可知W（ω），由方程（A.38）可知U（ω）。

问题全在于φ（ω）等于什么。

这个要求是当对W（ω）作傅氏反变换时，它将产生最小相位子波。

由方程（A.38）计算出函数U（ω），我们再对它作傅氏反变换回到时间域。

对它应用希尔伯特变换并使它返回到变换域，则：

W（ω）＝exp［UH（ω）］（A.42）

式中UH是U的希尔伯特变换。

这是保证子波最小相位的必要条件。

这个结论的精辟分析由Claerbut（1976）论述过了。

我们真正要知道的是子波的逆，因为这个子波的逆是脉冲反褶积的因子。

根据振幅谱和相位谱重新写（A.42）式，则有：

W（ω）=A（ω）ei

]（A.43）

这样，在傅氏变换域中的反褶积滤波是：

F（ω）=1/W（ω）=Af（ω）ei

]（A.44）

（A.43）式代人上式，有：

Af（ω）=1/Aw（ω）（A.45a）

φf（ω）=-φw（ω）（A.45b）

对（A.43）式作傅氏反变换得到的反褶积因子。

根据白噪化的反射系列的假设，我们可以把方程（A.37）中地震道的自相关函数代入。

我们通过图12流程图来概述频率域反褶积方法。

除了频率域处理方法在频谱计算中用于截断误差而产生微弱的噪音外，这两种方法实际上没有什么差别。

至此，由上述讨论可知，是等效的最小相位震源子波的逆。

通常所估算等效的最小相位震源子波的方法叫做谱的分解。

一旦计算出最小相位子波，那么对它求过得到反褶积因子就是一个简单事情了。

参照（A.45a），应该提醒，在振幅谱有“缺空”情况下就可能出现被零除的危险。

为保证滤波是稳定的，我们可以考虑在作除法之前把很小的数加到振幅谱上。

这就是预白噪化的基本概念。

那么方程（A.45a）的形式就变为：

Af（ω）=1/[Aw（ω）+ε（A.46）

此外，将振幅谱平滑以便使给出更稳定的因子。

另一个有效的方法是零相位反褶积。

这个反褶积可以通过选择方程（A.45b）给出的相位谱调整到零，这就有效地获得一个输入地震道记录展平振幅谱的反褶积算子而不改变它的相位。

这样的过程相当于时变谱白化，它起到在时间域中类似的作用。

A.5最佳维纳滤波

如下描述的是普通滤波模型。

我们希望设计一个滤波器f（t），使它的实际输出和期望输出之间的最小平方误差最小。

我们定该误差I为：

I＝

（A.47）

实际输出是滤波因子与输入的褶积；

yt＝ft*xt（A.48）

把方程（A.48）代入到方程（A.47）中，则：

（A.49）

这个目的是计算滤波因子时[f（0），f

（1），…，f（n-1）]使误差达到最小。

当然，滤波长度n预先确定。

数学上，我们规定I对f（i）的变化率为零：

=0i＝0，1，2，…，（n-l）（A.50）

将方程（A.49）平方项展开，有

I=

（A.51）

取偏导数为零：

=-2

（A.52）

或

i=0，1，…，（n-1）（A.53）

所以

（A.54）

最终，对每一个第i项，我们得到：

（A.55）

以矩阵形式表示为：

（A.56）

这里r是输入自相关的延迟，g是输入与期望输出之间的互相关。

我们可以应用莱文森递堆公式求出最佳维纳滤波因子f（i）。

现在我们计算包括在这个过程中所包含的误差，（A.51）式改写成下面形式：

Imin=

（A.57）

将下面的关系式

（A.58a）

（A.58b）

代人到（A.57）方程中，得到

（A.59）

（A.60）

最终

（A.61）

现在我们专门研究一下期望输出是输出时间超前型的模型d（t）＝X（t＋a）。

这样，互相关函数g变为：

（A.62）

应用定义

（A.63）

对于a＋τ，我们有

（A.64）

代人到方程（A.56）中，我们得到的求解预测滤波因子[f（0），f

（1），…，f（n-1）]的正规方程组。

（A.65）

现在考虑单位预测距离的特殊情况，方程（A．65）取下面的形式：

（A.66）

按等式右边，增广方阵的左边得到；

（A.67）

增加一行，并将滤波列上置入负号。

（A.68）

这样我们有n＋1个方程，n＋1个未知数，即[f（0），f

（1），…，f（n-1），V]。

实际上，根据方程（A.68），

V=r0-r1f0-r2f1-…-rnfn-1（A.69）

方程（A.68）与方程（A.36）是等效的，所以我们的结论是单位预测距离的预测滤波器和长度为（n＋1）的预测滤波器对于同样长度的反滤波是等效的。

应用方程（A.61），我们能够计算与单位预测距离滤波器有关的最小误差。

dt＝xt＋l，

（A.70）

并且还有

gt＝rt+1（A.71）

Imin=r0-（r1f0+r2f1+…+rnfn-1）（A.72）

它与方程（A.69）给出的数值V相同。

所以，当我们解方程（A.68）时，除了滤波系数之外，实际上是计算最小误差。

方程（A.65）给出了预测距离为a的预测滤波器。

这个期望输出是x（t＋a），我们说实际输出

（t＋a）是前者的估算值。

这个误差序列包括输入序列的未预测分量，定义为：

et＋a=xt+a-

（A.73）

我们假设未预测分量是我们希望从地震记录中提取反射率e（t）。

取方程（A.73）的z变换：

EaE（z）＝EaX（z）-F（z）X（z）（A.74）

E（z）＝[1-z-a-F（z）]X（z）（A.75）

现在我们确定一个新的滤波因子a（t），它的z变换是：

A（z）＝1-z-aF（z）（A.76）

所以，当代入到（A.75）式时，则有：

E（z）=a（z）X（z）（A.77）

和时间域的对应关系是

e（t）=a（t）*x（t）（A.78）

当确定e（t）后，则有：

反射率=a（t）*x（t）（A.79）

方程（A.79）说明的是如果我们把滤波因子a（t）应用于输入地震道，我们就可以获得反射率序列，这恰好是我们做反褶积的目的。

这个滤波因子a（t）是由方程（A.76）给的，在时间域里有：

at=（

）（A.80）

我们根据预测滤波因子f（t）得到滤波因子则a（t），f（t）是方程（A.65）的解。

我们称a（t）为“预测误差滤波器”，在单位预测距离的特殊情况下，由方程（A.80）给出的滤波因子取下列形式：

at＝（l，-f0，-f1,-f2，…，-fn-1）（A.81）

又可看出，这与方程（A.68）的解是一样的。