反射系数反演原理

反射系数反演原理提高地震资料分辨率一直是地震资料数字处理工作者的追求目标。

高分辨率地震资料数据处理的关键环节就是压缩地震子波，或者去除地震波在地下传播过程中干涉、调谐等效应对地下地层的影响，拓宽有效地震信号的频带范围，特别是较为准确地拓宽高频成分。

薄层的识别和厚度估计是当前地震石油勘探领域的主要研究方向。

薄层反射系数反演方法也就应运而生。

早在1999年Partyka在发表的文章中提出应用时频分析方法计算薄层厚度，即为薄层陷频法，他认为少数几个反射系数谱不像是长时窗反射系数谱一样是白噪，在频率域存在周期性的陷频规律，地震道的谱就是反射系数谱与子波谱的乘积，消去子波谱就能显示出陷频谱，从而进行薄层厚度的计算。

2001年K.J.Marfurt提出用滑动时窗的频谱分析方法计算多种频率相关的属性进行储层厚度定性解释，对薄窄河道展布规律进行描述。

上述两种方法尝试利用频谱分解的结果来求取薄层厚度，但如果地震频带宽度不足以清晰识别谱峰和陷频变化规律时，谱分解对于分辩薄层还是存在困难，这也推动了新方法的发展，无需精确识别频宽内的波峰和波谷，即为谱反演方法。

2005年Portniaguine和Castagna提出了一种叠后谱反演方法，来解决在小于调谐厚度时的薄层预测问题。

这个方法更多的从地质上去考虑，而不是数学上的假设。

其重点在于通过分频方法来获取局部频谱信息。

这种谱反演或薄层系数反演方法最终输出的为反射系数序列，其视分辨率要远高于输入的地震数据，可以用来对薄储层进行精细的描述和刻画。

并指出了该方法具有不需要任何先验模型、反射系数的数学假设、层位约束，也不需要井资料强制约束等优点，并可用来分辨小于调谐厚度的薄层[5,6]。

之后很多学者在谱反演方面做出很多研究成果。

Puryear和Castagna (2008)对谱反演理论给出了详细的说明，主要包括Widess楔形模型理论，把反射系数序列分解成偶分量和奇分量，发展了一种新的谱反演算法[7]。

Satinder Chopra, John Castagna和Yong Xu(2009)由谱反演所得的反射系数算出了波阻抗剖面，并将其应用到了层厚的确定和地层学解释中[8,9]。

袁三一博士和王尚旭教授等（2009）提出了一种相对快速的谱反演系数反演混合技术，其采用求解的算法是Particle Swarm Optimization (PSO)，粒子群算法和Levenberg - Marquardt(L-M)[10,11]。

Kelyn Paola Castaño 和 Germ án Ojeda(2010)采用遗传算法和模拟退火优化算法来求解Castagna 教授的基于反射系数序列奇偶原理的谱反演目标函数，并指出遗传算法效果比模拟退火算法稍好[12]。

1. Widess 模型1973年的薄层反射Widess 模型认为地震分辨率极限为1/8波长，如图1中所示，Widess 模型与30Hz 雷克子波褶积后峰值频率与层厚关系。

频率域内地震频谱振幅在1/4波长达到最大，当层薄至1/8波长时地震波形及频率不会再有明显变化，振幅逐渐减小。

实际情况下，假设条件是难以保证，顶底反射系数比不可能是连续变化的，不同反射系数比下峰值振幅、峰值频率与厚度的关系也不尽相同如图2。

图3奇偶分量峰值振幅、峰值频率与厚度的关系我们利用Chung和Lawton(1995)提出的方程计算了峰值振幅和峰值频率，图3分别表示为随层厚减薄时由奇部、偶部形成的反射系数对模型峰值振幅及峰值频率的变化。

根据这种模型的分析，显示了在偶部序列占优时奇部序列和偶部序列主峰振幅和频率的随层厚的变化特征，同时也显示了总的主频变化特征。

当薄层的顶底界面处反射系数并不相等、极性也不相反时，可以看到在小于调谐频率的一半时主频会随着层薄而减小。

这种变化的转折点位置取决于序列中偶部和奇部反射系数的相对大小。

与Widess模型相反的是，在层厚小于该点后，主频还依赖于地层厚度。

这说明与以往观点对比，地层响应对薄层更敏感。

由于奇偶分量峰值振幅、峰值频率随层厚的变化规律的确定性，地震资料的主频确定的情况下，对于单层厚度的地震模型，可以通过奇偶分解的方式，准确预测小于调谐厚度下的薄层厚度。

反射系数经过奇偶分解之后，奇分量与偶分量的主频随地层厚度变化是一个确定关系，通过这种关系，在主频已知的情况下可进行厚度的预测。

图4 单层厚度模型奇偶分解预测厚度图3. 谱反演的目标函数谱反演的原理就是根据时间域褶积模型，从地震记录中去除地震子波的影响，进而得到反射系数序列。

在时间域内一个脉冲对的表达为(如图5)：1121()()()g t r t t r t t T δδ=-+-- (1)其中1r 为层顶部反射系数，2r 为层底部反射系数，t 为时间位置，1t 为顶部反射的时间位置，T 为层厚度。

图5 两层反射率模型将分析点放在层的中心点位置，进行傅立叶变换后，用三角法则进行化简约去t ，得到：()2cos()2sin()e o g f r fT r fT ππ=+ (2)基于褶积原理，),(f t s 是地震数据，),(f t w 是已知的子波。

()⎰-+=wwt t o e dt t fT t ir t fT t r f t w f t s ]}sin[)()](cos[)({),(),(ππ (3)最后整理得到单层目标函数：()()()()()()()(),,,Re ,/,cos Im ,/(,)sin wHw wl w e o t e e f t t f o o t O r r T t a S t f W t f r t fT t dt df a S t f W t f r t fT t dt ππ-=⎡⎤⎧⎫⎪⎪⎢⎥-+⎡⎤⎡⎤⎨⎬⎣⎦⎣⎦⎢⎥⎪⎪⎩⎭⎢⎥⎧⎫⎢⎥⎪⎪-⎡⎤⎡⎤⎨⎬⎢⎥⎣⎦⎣⎦⎪⎪⎢⎥⎩⎭⎣⎦⎰⎰⎰ (4) 目标函数中，()(),/,S t f W t f 可以认为是反褶积过程，在这里我们假设子波是已知。

子波的正确与否将很大程度上影响反演的结果，因此子波的提取是一个非常重要的过程。

但是我们对于子波往往是不知道的，如果子波知道了，反射系数也就知道了。

因此我们做的是尽可能的使子波接近于真实子波。

厚度T,我们在反射系数中可以认为是第一层和最后一层的距离，第二层和倒数第二层的距离，依次类推。

在信号的奇偶分解中 我们可以认为是进行奇偶分解的信号对之间的距离。

以上的推导都是基于窗口中心位置的坐标为零，但实际上数据都是从零开始的，而窗口的中心位置坐标不可能为零，而是半个窗口的位置。

因此我们实际应用中应该对上面的推导加一时移。

实部和虚部时移后的谱可以表示为：2Im[()]2sin()cos(2)2cos()sin(2)i f t o e e g f r fT f t r fT f t πππππ∆=∆+∆ (5) 2Re[()]2cos()cos(2)2sin()sin(2)i f t e o e g f r fT f t r fT f t πππππ∆=∆-∆ (6)多层的目标函数变为：()()()()()(){}()()()(){}002()cos cos22sin 22sin cos22()c ,,,Re ,/,sin Im ,/(,)sin os 2wHwl e o e o t e f t f o o e O r r T t a S t f W t f fT t f t fT t f t dt df a S t f W t f fT t f t fT t r t r t r t r t f t dt ππππππππ+⎡⎤⎡⎤⎣⎦⎣⎦⎡⎤⎡⎤=⎡⎤⎧⎫⎪⎪-∆∆+⎡⎤⎢⎥⎨⎬⎣⎦⎪⎪⎢⎥⎩⎭⎢⎥⎧⎫⎪⎪⎢⎥-∆-∆⎡⎤⎨⎬⎣⎦⎢⎥⎪⎪⎩⎭⎣⎦⎦⎣⎣⎦⎰⎰⎰(7)其中，w t 为半个窗口的宽度；e a ，o a 为平衡参数。

其矩阵形式，形如：1111212233112233(()/())cos(2t )cos(2t )cos (()/())(()/())(()/())(()/())(()/())(()/())(()/())m m m m real X f W f f f real X f W f real X f W f real X f W f imag X f W f imag X f W f imag X f W f imag X f W f ππ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥-⎢⎥⎢⎥-⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦31n 1122232n 2132333n 3123n 112131(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )cos(2t )sin(2t )sin(2t )sin(2t )sin(2m m m m f f f f f f f f f f f f f f f f f ππππππππππππππππππ123n 1122232n 2n 132333n 3123n (t )(t )(t )t )sin(2t )sin(2t )sin(2t )sin(2t )(t )sin(2t )sin(2t )sin(2t )sin(2t )sin(2t )sin(2t )sin(2t )sin(2t )m m m m r r r f f f f f r f f f f f f f f ππππππππππππ⎡⎤⎢⎥⎢⎥⎢⎥⎡⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦(8)通常地震道采样点数远大于反射系数的个数，因此谱反演理论一部分就是为了解决(8)式的超定线性方程组求解的适定性和稳定性问题。

通常的迭代计算流程如下：图6 反射系数迭代计算流程任何反演方法都是一个系统工程，许多小环节的处理手段都会影响到反演的结果，子波提取、噪声及初始模型选取都对反演结果有影响。

在实际地震数据的处理中，反演过程大致包含以下几步：步骤一、子波提取，精细提取时变子波；步骤二、从地震数据中去除子波，提取计算反射系数的奇部和偶部；步骤三、根据稀疏脉冲反演计算初始反射系数的方法，反演出高频成分；步骤四、高频成分与奇部、偶部反射系数，由权重函数控制，组合出完整的宽频反射系数体。