2017, Vol. 60
2. 中海石油(中国)有限公司深圳分公司, 深圳 518000;
3. 中国地质大学(北京), 北京 100083;
4. 北京京能油气资源开发有限公司, 北京 100022
2. CNOOC Ltd. _Shenzhen, Shenzhen 518000, China;
3. China University of Geosciences(Beijing), Beijing 100083, China;
4. Beijing Energy Oil & Gas Resources Development Co. Ltd., Beijing 100022, China
近些年来，非常规油气资源的勘探开发在业界备受关注，其油气产量占油气总产量的比例快速升至10%以上.不同于以圈闭评价为主的常规油气资源，非常规油气资源勘探的目标是突破“甜点区”，确立连续性或准连续性油气边界.而地层的脆性指数，作为“甜点区”评价标准的六特性之一，对预测甜点范围、指导水力压裂均具有重要的意义.
脆性指数的计算既可以通过脆性矿物(如石英)的体积分数来表征，也可以借助于杨氏模量、泊松比等弹性参数组合表示.因此可以利用常规叠前地震反演获取的P波阻抗(速度)、S波阻抗(速度)和密度，将其转换为杨氏模量和泊松比，并进行归一化，进而得到脆性指数.考虑到间接计算带来的累计误差，宗兆云等(2012)对Aki-Richards近似式进行重新推导，得到关于杨氏模量、泊松比和密度的YPD近似式，在贝叶斯叠前反演框架下，可以直接提取稳定的杨氏模量和泊松比.
考虑到基于AVO近似式的线性反演具有诸多假设——中小入射角、常背景纵横波速度比以及弱弹性参数反射率等，本文从精确Zoeppritz方程出发，将其直接表达为关于脆性指数、P波速度和S波速度的函数(BI_Zoeppritz方程)，借助于广义线性反演技术可以直接获取脆性指数.由于需要同时反演上下地层的脆性指数、P波速度和S波速度，Jacobian矩阵的条件数较大，为了降低反演的不适定性，我们采取“两步法反演”，第一步对不同入射角的部分叠加剖面采用基追踪-稀疏层反演技术获取高分辨率角度反射系数，消除角度子波对后续反演的影响；第二步对已获取的角度反射系数进行BI_Zoeppritz方程广义线性反演，便可以直接提取高分辨率、高精度的脆性指数.
2 基于脆性指数的精确Zoeppritz方程表达Zoeppritz方程作为AVO技术及叠前地震反演的核心，在已知上下地层的P波速度、S波速度和密度情况下，可以计算不同入射角的反射系数和透射系数.张丰麒等(2013)、Fang等(2016)为了减少反演参数的个数，将Zoeppritz方程表达为关于P波反射率、S波反射率、密度反射率以及上层纵横波速度比这四个参数的函数，但是反演中需要假设上层纵横波速度比为背景纵横波速度比；魏修成等(2012)将Zoeppritz方程表达为关于射线参数和弹性参数的函数，并提出了基于射线参数道集的AVO反演方法.为了直接建立脆性指数与角度反射系数之间的关系，利用弹性参数之间的转换关系，可以将精确Zoeppritz方程表达为上下地层脆性指数、P波速度和S波速度的函数.考虑到脆性指数的表达式为：
|
(1) |
(1) 式中的BIYM和BIPR分别表示归一化处理后以百分数表示的岩石杨氏模量和泊松比，其取值范围为[0, 100]，这主要是考虑到杨氏模量和泊松比的量纲差异较大，为了均衡反映它们在评价岩石脆性中的作用而做出的数据处理.其表达式如下：
|
(2) |
(2) 式中YM和PR分别表示储层的杨氏模量和泊松比，YMmin和YMmax分别表示目的层段最小和最大杨氏模量，PRmin和PRmax分别表示目的层段最小和最大泊松比.结合杨氏模量、泊松比与P波速度、S波速度和密度的关系：
|
(3) |
可以直接推导出脆性指数和P波速度、S波速度以及密度的关系：
|
(4) |
由于YMmin, YMmax, PRmin, PRmax可以从实际工区的先验信息中获取，是已知参数.因此脆性指数BI可以看成P波速度和S波速度的非线性函数，密度的线性函数.对(4) 式进行整理，可以得到密度关于脆性指数、P波速度和S波速度的函数：
|
(5) |
为了后续表达的方便，将公式(5) 简写为：
|
(6) |
将精确Zoeppritz方程中的密度项全部替换为(5) 式：
|
(7) |
(7) 式称之为BI_Zoeppritz方程.式中α1, α2, β1, β2分别表示PP波反射角、PP波透射角、PS波反射角和PS透射角，通过BI_Zoeppritz方程可以直接建立地层脆性指数与PP波反射系数之间关系.
3 脆性指数反演的可行性分析为了检验BI_Zoeppritz方程的精度，通过以下双层模型计算反射系数随角度的变换关系，并与精确Zoeppritz方程进行对比.其模型如表 1所示.
|
表 1 双层模型 Table 1 Two-layer model |
假设最小杨氏模量YMmin=7×109、最大杨氏模量YMmax=6×1010、最小泊松比PRmin=0.1和最大泊松比PRmax=0.4.根据公式(4) 可以计算出层一的脆性指数为33.5836，层二的脆性指数为34.7656，然后分别利用BI_Zoeppritz方程和精确Zoeppritz方程计算1°到90°的PP波反射系数和PS波反射系数.其计算结果如图 1所示，可以看出无论是PP波反射系数还是PS波反射系数，在1°到90°范围内BI_Zoeppritz方程与精确Zoeppritz方程的计算的反射系数曲线是完全重合的，这是由于BI_Zoeppritz方程的推导过程并未引入任何假设，因此理论上是没有精度损失的.
|
图 1 精确Zoeppritz方程与BI_Zoeppritz方程精度对比 (a) PP波反射系数；(b) PS波反射系数. Fig. 1 Comparison of exact Zoeppritz equation and BI_Zoeppritz equation (a) PP wave reflectivity; (b) PS wave reflectivity. |
为了进一步检验BI_Zoeppritz方程中PP波反射系数对脆性指数、P波速度与S波速度的敏感性，可以借助反射系数对各项参数的一阶偏导数来表征不同参数对反射系数的贡献率.在相同角度的情况下，一阶偏导数越大，说明反射系数对该参数的变化越敏感，该参数对反射系数的贡献量较大，不确定性较小.
考虑到BI_Zoeppritz方程中反射系数是关于上、下两层介质参数的函数，并且脆性指数与P波速度、S波速度的量纲不一致，因此通过计算反射系数对各参数的变化率的一阶偏导数代替更为合理：
|
(8) |
(8) 式中
|
图 2 弹性参数模型(黑线)以及对应的初始模型(红线) Fig. 2 Elastic model (black line) and corresponding initial model (red line) |
图 3给出了入射角分别为5°、20°以及35°时，反射系数对脆性指数变化率(蓝线)、P波速度变化率(红线)以及S波速度变化率(绿线)示意图.可以看出在小角度的情况下，反射系数对P波速度变化率的一阶偏导数最大，S波速度其次，脆性指数最低，说明在小角度情况下，P波速度对反射系数的贡献率最大；随着角度增加到20°，反射系数对脆性指数变化率的一阶偏导数逐步增大，基本与S波速度持平；当角度为35°的情况下，反射系数对脆性指数变化率的一阶导数继续增加，基本与P波速度持平.可以看出，随着角度的增大，反射系数对脆性指数变化的敏感程度逐步提高，因此对反射系数的贡献率也是逐步增加.那么借助BI_Zoeppritz方程来反演脆性指数，需要有足够大的入射角，才能保证脆性指数反演的精度和稳定性.
|
图 3 反射系数对各参数变化率的一阶偏导数 (a)入射角为5°；(b)入射角为20°；(c)入射角为35°. Fig. 3 First-order derivatives of ratio of reflectivity to each parameter (a) Incident angle 5 degrees; (b) Incident angle 20 degrees; (c) Incident angle 35 degrees. |
从BI_Zoeppritz方程可以看出，PP波反射系数RPP是关于BI1, BI2, VP1, VP2, VS1, VS2这六个参数的非线性函数，即
|
(9) |
对于非线性函数，可以利用全局寻优算法求解，但是计算量会非常大.另外还可以通过广义线性反演的手段，在给定初始模型的情况下进行迭代求解.李录明等(2010)利用广义线性反演成功地从叠前地震数据中获取了弹性参数和各向异性参数；Huang等(2015)将贝叶斯理论和广义线性反演相结合，通过引入合适的正则化项，可以有效降低广义线性反演的不适定性；张丰麒等(2013)、Fang等(2016)结合三变量柯西分布和广义线性反演分别实现了基于精确Zoeppritz方程的常规P波反演和多波联合反演.
在初始模型已知时，(9) 式可以利用泰勒一阶展开式进行线性展开：
|
(10) |
(10) 式中RPP0表示由初始模型计算获取的PP波反射系数，尽管需要求解的参数多达六个，但是考虑到第i个地层分界面对应的下层弹性参数为第i+1个地层分界面对应的上层弹性参数，同时第i个地层分界面对应的上层弹性参数又是第i-1个地层分界面对应的下层弹性参数.因此当扩展到N个地层时，由于上下层是相互耦合的，仅需要三个角度便可以求解出来，其正演矩阵表达式为：
|
(11) |
(11) 式中Gx(θi)表示第θi个入射角下，由
|
(12) |
ΔRPP(θi)表示观测的第θi角度反射系数向量与初始模型计算的角度反射系数向量的差.因此借助(11) 式进行迭代计算便可以获取高精度的反演参数.但是当观测数据中包含噪声时，由于Jacobian矩阵的条件数较大，反演异常不稳定，通过借助贝叶斯理论引入合适的正则化项，可以提高广义线性反演的稳定性.考虑到ΔBI, ΔVP, ΔVS的先验分布无法从测井中获取，可以将(11) 式重组为反演参数为BI, VP, VS的矩阵表达式(张丰麒等，2013；Fang et al., 2016)：
|
(13) |
其中：
|
(14) |
为了后续表达方便，将(13) 式简写为：
|
(15) |
借助于贝叶斯定理，假设角度反射系数中的噪声服从多变量高斯分布，同时模型参数BI, VP, VS也服从多变量高斯分布，利用贝叶斯定理整合角度反射系数的似然函数和模型参数的先验分布得到：
|
(16) |
(16) 式中Cn为噪声的协方差矩阵，考虑到不同道的噪声相互独立，因此Cn=σn2I，σn2为噪声的方差，I为单位阵，Cx为模型参数的协方差矩阵，可以利用公式(17) 从测井数据中获取：
|
(17) |
(17) 式中cx(i, j)表示3×3矩阵cx的第i行，第j列元素，其计算公式如下：
|
(18) |
xtrend表示模型参数的低频趋势，可以通过对测井曲线进行低通滤波获取.
利用模型参数后验概率极大等价于反演目标函数极小值，将(16) 式表达为反演目标函数：
|
(19) |
(19) 式中的μ为超参数，在控制反演精度和先验约束权重上起到折中的作用.对(19) 式两边求导，并令其偏导数为0，推导出：
|
(20) |
当给定初始模型x0，首先计算Jacobian矩阵G，然后利用(13) 式计算H，接着利用(20) 式计算首次迭代的模型参数x，将获取的模型参数x作为下次反演的初始模型，重新计算新的模型参数，经过3~5次迭代便可以获得高精度的反演结果.
值得注意的是，借助于(20) 式获取脆性指数的前提是观测的角度反射系数RPP是已知的，并且Jacobian矩阵G中并未包含子波项，这主要是考虑到如果在G矩阵中融合子波褶积矩阵，G矩阵的不适定性将更强(叠前地震反演的不适定性有两方面因素：1) 子波的带限特性导致子波矩阵的逆是非唯一，这是地震反演问题非唯一的根本；2) 由弹性参数向角度反射系数的转换的正演矩阵的条件数较大，在含有噪声的情况下，其解也是不稳定的)，并且无法获取高分辨率反演结果.
综合上述分析，我们采用分步反演思想，首先利用基追踪-稀疏层反演将角度部分叠加数据转化为角度反射系数，然后再采用BI_Zoeppritz方程广义线性反演提取脆性指数.
5 基追踪-稀疏层反演提高地震垂向分辨率一直是当前地震勘探领域内的一个重点和难点.Widess(1973)提出地震的极限分辨率为1/8波长，在考虑噪声的时候，实际上只能分辨到1/4波长.Castagna(2004)，Chopra等(2009)基于楔形模型提出反射系数奇偶分解理论——一个反射系数对可以唯一分解为一个奇分量和一个偶分量，并认为奇分量不利于薄层的识别，而偶分量则有助于薄层的识别；在此基础上，Partyka(2005)、Portniaguine和Castagna(2005)、Puryear和Castagna(2008)、Chopra等(2009)借助谱分解和模拟退火实现了高分辨率的谱反演；Nguyen和Castagna(2010)通过最小化模型参数的L0范数，利用匹配追踪将地震道分解为奇分量和偶分量系数；Chen等(2001)指出基追踪相对匹配追踪计算效率更高，并且当过完备字典的原子不是正交的情况下，依然相对稳定；Zhang和Castagna(2011)随后提出基追踪-稀疏层反演算法，在此基础上发展了基于空间约束的多道基追踪反演(Zhang et al., 2013a)以及叠前AVA基追踪反演(Zhang et al., 2013b)；刘晓晶等(2016)、印兴耀等(2016)发展了基于基追踪的弹性阻抗反演和基于模型约束的基追踪反演.
Castagna(2004)和Chopra等(2009)提出一个反射系数对可以唯一分解为一个反射系数奇分量和一个反射系数偶分量，如图 4所示(Zhang and Castagna, 2011).
|
图 4 反射系数奇偶分解示意图(Zhang and Castagna, 2011) Fig. 4 Schematic diagram of dipole decomposition of reflectivity (Zhang and Castagna, 2011) |
图 4中r1, r2构成一个反射系数对，反射系数奇分量由符号相反的两个狄拉克δ函数构成，而偶分量则由符号相同的两个狄拉克δ函数构成；a，b分别表示偶分量和奇分量的系数，利用如下数学表达式可以表达反射系数奇偶分解的过程：
|
(21) |
(21) 式中re和ro分别表示反射系数的偶分量和奇分量，N表示反射系数序列的采样点数，M表示一个薄层可能的最大时间采样点数，Δt代表时间采样率，则MΔt表示一个薄层可能的最大时间厚度.(21) 式可以表达为矩阵形式：
|
(22) |
其中的m表示奇偶分量对应的系数：
|
(23) |
而D表示反射系数奇偶分解矩阵，以采样点N=4为例，D的表达式如下：
|
(24) |
|
(25) |
(24)、(25) 式的地球物理含义是，当地层有4个采样点时，薄层可能的厚度为Δt、2Δt和3Δt，其中当i=1时，由3组反射系数对构成，意味着厚度为Δt的薄层在地层中位置会有3种情况；当i=2时，由2组反射系数对构成，意味着厚度为2Δt的薄层在地层中位置会有2种情况；当i=3时，由1组反射系数对构成，意味着厚度为3Δt的薄层在地层中位置只有1种情况.据此推断出，当地层时间采样点为N时，最多可达(N－1)+(N－2)+…+(N－M)组反射系数对，每一组反射系数对则对应一个可能的薄层.因此反射系数奇偶分解矩阵D包含了所有的薄层厚度和薄层位置的信息，该信息是冗余过完备的.然而在针对实际数据进行处理时，薄层最大的时间厚度MΔt不宜过小或者过大.当MΔt过小时，意味着包含的薄层信息不够充分，无法有效恢复真实的反射系数；当MΔt过大时，会造成D矩阵维数急剧增加，因此计算效率显著降低.通常MΔt应该选择1/4波长附近，即
结合地震褶积理论和(22) 式，则地震正演过程表达如下：
|
(26) |
(26) 式中，W为子波褶积矩阵，s为实际地震信号，n表示随机噪声.在模型参数L1范数极小化的约束下，稀疏层反演的目标函数为：
|
(27) |
(27) 式正好对应于信号稀疏表示领域中的基追踪问题，基追踪可以在楔形字典WD中寻找最优正交基，而地震信号s在这个正交基上的投影便对应于模型参数m的非零元素，这些非零元素对应于薄层顶、底的奇偶分量系数.常规地震反演以反射系数为研究对象，而基追踪-稀疏层反演以薄层的厚度和空间位置为研究对象.
针对基追踪问题的求解，由于奇偶分量系数m≥0，可以利用不等式约束的最优化方法中的内点法或可行方向法进行求解.其中内点法将不等式约束问题转换为序列最小二乘优化问题，例如原始-对偶对数障碍法等；而可行方向法是指寻找一个使目标函数值减小的同时又使得自变量始终在有效集范围内的可行性方向，例如梯度投影稀疏重构等.
Zhang等(2013a)指出利用(27) 式直接进行反演，其横向连续性较差.为此，提出了基于多道的稀疏层反演.然而多道反演无疑需要消耗更大的内存，为此本文提出在(27) 式中引入低频软约束项，测井低频信息的融入可以有效增强反演剖面的横向连续性.
|
(28) |
(28) 式中第二项低频软约束可以这样理解：由于Dm=r，C为积分矩阵，则CDm等价于对角度反射系数沿垂向积分，因此等价于弹性阻抗的自然对数；L表示低通滤波矩阵，主要由正反离散傅里叶变换矩阵和基于汉明窗的低通滤波对角阵构成，具体表达式见附录，因此LCDm表示反演获取的弹性阻抗的自然对数的低频成分，EItrend表示由井曲线计算获取的弹性阻抗的自然对数的低频趋势，通过该项约束可以保证两者在L2范数意义下差值最小，使得反演获取的低频信息同测井低频信息相吻合.
将基追踪-稀疏层反演同基于精确Zoeppritz方程脆性指数反演相结合，图 5给出了基于基追踪-广义线性脆性指数直接反演的流程.
|
图 5 基于基追踪-广义线性脆性指数直接反演流程 Fig. 5 Flowchart of direct inversion for brittleness index based on GLI with basic-pursuit decomposition |
合成角道集的建立是通过图 6中的测井曲线，在不考虑地震波透射损失、几何扩散以及多次波的情况下，利用精确Zoeppritz方程分别计算15°、30°和45°的角度反射系数，然后再和30 Hz主频的雷克子波褶积得到，如图 7a所示，图 7b给出了信噪比为4的含有随机噪声的合成角道集.
|
图 6 弹性参数模型及其初始模型(黑线为井曲线，绿线为初始模型) Fig. 6 Elastic parameter model (black) and initial model (green) |
|
图 7 (a)合成角道集；(b)含噪声合成角道集 Fig. 7 (a) Synthetic angle gather; (b) Synthetic angle gather with noise |
为了展示基追踪-BI_Zoeppritz方程广义线性脆性指数反演算法的优势，分别进行了如下三种反演：常规叠前AVO反演、基于BI_Zoeppritz方程的广义线性同步反演(指在(20) 式的G矩阵中融入子波褶积矩阵后，直接利用(20) 式进行反演)以及基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演.其中常规叠前AVO反演的三参数为P波速度、S波速度以及密度，因此脆性指数需要经过间接转换才能获得.基于BI_Zoeppritz方程广义线性同步反演和基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演的初始模型是通过对测井曲线进行低通滤波获取(如图 6中的绿线所示).图 8、图 9和图 10分别给出了这三种反演结果与井曲线的叠合显示.
|
图 8 常规叠前AVO反演结果 (黑线为实测井曲线，红线为反演结果) Fig. 8 Real log curve (black) and routine pre-stack AVO inversion (red) |
|
图 9 基于BI_Zoeppritz方程广义线性同步反演结果 (黑线为实测井曲线，红线为反演结果) Fig. 9 Real log curve (black) and inversion result of GLI based on BI_Zoeppritz equation (red) |
|
图 10 基于基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演(黑线为实测井曲线，红线为反演结果) Fig. 10 Real log curve (black) and direct inversion of brittleness index based on GLI with basic-pursuit decomposition (red) |
由于常规叠前AVO反演是基于稀疏脉冲反演算法，其垂向分辨率较低，并且脆性指数的反演结果与井曲线结果吻合度较差(图 8).相比常规叠前AVO反演算法，基于BI_Zoeppritz方程广义线性同步反演，由于采用了精确Zoeppritz方程并且基于迭代算法，因此脆性指数的反演精度得到了提高，但是由于子波的带限特性对地震反演的影响，其垂向分辨率仍然较低(图 9).基追踪-稀疏层反演综合考虑了薄层的厚度和位置信息，可以有效提高反演的垂向分辨率.在分别获取高分辨率的三个角度反射系数(图 11)后，利用BI_Zoeppritz方程广义线性反演对角度反射系数进行分解，进而可以获得高分辨率、高精度的脆性指数(图 10).为了进一步量化三种反演获取的脆性指数的精度，分别统计了其脆性指数反演结果与实际脆性曲线的相干系数，如表 2所示.
|
图 11 利用基追踪-稀疏层反演获取的三个角度反射系数(黑线为实测井曲线，红线为反演结果) Fig. 11 Real log curve (black) and three reflectivities from sparse-layer inversion based on basic-pursuit decomposition (red) |
|
|
表 2 三种反演获取的脆性指数与实测脆性曲线相干系数统计 Table 2 Statistics of correlation coefficients for three inverted brittle indexes and real log curves |
在模型试算过程中，本文分别对基于BI_Zoeppritz方程广义线性同步反演和基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演的Jacobian矩阵的条件数进行了统计，如表 3所示.
|
|
表 3 不同反演对应的Jacobian矩阵条件数统计 Table 3 Statistics of conditional numbers of Jacobian matrix in two inversions |
从表 3中可以看出基于BI_Zoeppritz方程广义线性同步反演的Jacobian矩阵条件数要明显高于基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演的，这主要是由于同步反演法中Jacobian矩阵耦合了子波褶积矩阵，子波的带限特性提高了反演的不适定性.
为了检验基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演的抗噪性，对图 7b中含噪声合成角道集进行反演.图 12给出了含噪声情况下的反演结果，相比无噪声的反演结果，反演精度有所下降，弱反射区域反演结果与实测曲线拟合度较差，主要是因为弱反射被随机噪声淹没，但从整体来看，反演精度尚可，可见基追踪-BI_Zoeppritz方程广义线性脆性指数直接反演具有较强的抗噪性.
|
图 12 含噪声合成道集的反演结果 (黑线为实测井曲线，红线为反演结果) Fig. 12 Real log curve (black) and inversion result of angle gather with the noise (red) |
实际资料取自A工区，该工区致密砂岩较为发育，为了对该工区进行脆性评价，需要借助岩石物理与叠前反演技术来预测脆性指数.在经过去除线性干扰、压制多次波、叠加速度分析后，通过真振幅叠前时间偏移后获取的共反射点(CRP)道集，并转换为角道集，再进行分角度部分叠加，其小角度、中角度和大角度部分叠加剖面分别如图 13a、图 14a和图 15a所示.
|
图 13 (a)小角度叠加剖面；(b)基追踪反演获取的小角度反射系数剖面 Fig. 13 (a) Stacked section of small angle; (b) Inverted reflectivity section of small angle |
|
图 14 (a)中角度叠加剖面；(b)基追踪反演获取的中角度反射系数剖面 Fig. 14 (a) Stacked section of medium angle; (b) Inverted reflectivity section of medium angle |
|
图 15 (a)大角度叠加剖面；(b)基追踪反演获取的大角度反射系数剖面 Fig. 15 (a) Stacked section of large angle; (b) Inverted reflectivity section of large angle |
在此基础上，对这三个角度部分叠加数据体进行井震标定，在统一时深关系下，分别提取了小角度、中角度和大角度子波.通过测井数据计算获取三个角度的弹性阻抗，并采用其0~10 Hz的低频信息进行约束，然后利用基追踪-稀疏层反演获取高分辨率的角度反射系数，如图 13b，图 14b和图 15b所示.从原始地震剖面和角度反射系数剖面对比可以看出，基追踪-稀疏层反演可以有效消除角度子波的旁瓣效应，使地层分界面更加清晰，提高了地震数据的垂向分辨率，并且角度反射系数剖面的信噪比和横向连续性也较高.
在获取高分辨率的角度反射系数剖面后，利用测井数据对目标层段的杨氏模量和泊松比的分布范围进行了统计，得出杨氏模量的最小值和最大值分别为7.84661×109和4.2913×1010，泊松比的最小值和最大值分别为0.0839925和0.414711，从而可以根据P波速度、S波速度和密度计算出对应的脆性指数曲线，进而为反演建立有效的先验约束.图 16a是利用常规叠前AVO反演通过间接转换获取的脆性指数反演剖面，图 16b给出了基追踪-BI_Zoeppritz方程广义线性脆性指数反演剖面.通过这两张图的对比可以看出，传统方法获取的脆性指数剖面受限于地震频段，其垂向分辨率和地震一致.而基追踪-BI_Zoeppritz方程广义线性脆性指数反演剖面垂向分辨率较高，展示了丰富的细节信息，并且在横向变化上也更加清晰明显.为了进一步观察传统方法与新方法的反演精度和垂向分辨率，图 17给出了井B(位于第76CDP附近)的井旁道反演结果与实测井曲线叠合对比.可以看出，新方法相对传统方法的三参数反演结果的垂向分辨率要高，且与井曲线的吻合程度较好，特别是在第1600~1700 ms附近.
|
图 16 (a)常规叠前AVO反演间接计算获取的脆性反演结果; (b)基追踪-BI_Zoeppritz方程广义线性脆性反演结果 Fig. 16 (a) Brittleness index section from routine pre-stack AVO inversion; (b) Brittleness index section from the direct inversion based on GLI with basic-pursuit decomposition |
|
图 17 井旁道反演结果与实测井曲线对比叠合显示 (黑线为井曲线；蓝线为传统叠前AVO同步反演结果；红线为基追踪-BI_Zoeppritz方程广义线性反演结果) Fig. 17 Comparison of inverted results from routine pre-stack AVO (blue) and GLI with basic-pursuit decomposition (red) and the log curve (black) |
本文根据弹性参数之间的转换关系，对Zoeppritz方程进行了重组，将其表达为脆性指数、P波速度和S波速度的函数——BI_Zoeppritz方程，从而建立起角度反射系数与脆性指数的直接关系.
考虑到常规广义线性反演需要在Jacobian矩阵中耦合子波褶积矩阵，造成反演的不适定性进一步增强，并且无法获取高分辨率的反演结果.本文通过利用基追踪-稀疏层反演对角度叠加数据进行处理，首先获取高分辨率的角度反射系数，然后再用广义线性反演对角度反射系数进行分解，可以直接提取高精度、高分辨率的脆性指数，通过模型试算和实际数据测试，验证了基追踪-BI_Zoeppritz方程脆性指数直接反演方法的有效性.
由于该反演方法需要针对每个角度道单独进行基追踪-稀疏层反演，因此不同角度道的反射稀疏约束项权重λ选择的不同，是否会影响不同角度之间反射振幅的相对关系仍需进一步研究.另外也可以针对不同角度开展基追踪多道同时反演，进一步保障不同角度之间反射振幅的相对关系.
附录
低频软约束中的低通滤波矩阵L的具体表达式：
|
(A1) |
(A1) 式中E矩阵起到将模型参数进行三倍延长的作用.考虑到反演目的层段时窗过短会引起在频域的分辨率不够，或者由于数据两端的截断效应会引起低通滤波结果不准确，通过该矩阵可以使得在数据两端各以原数据的倒序进行延长，E矩阵的数学表达式为：
|
(A2) |
F和F-1分别为DFT正反变换矩阵，起到将数据进行傅里叶正、反变换的作用.
|
(A3) |
|
(A4) |
(A3) 式和(A4) 式中的
Castagna J P. 2004. Spectral decomposition and high resolution reflectivity inversion.//Presented at the Oklahoma Section Meeting. SEG.
|
|
Chen S S, Donoho D L, Saunders M A.
2001. Atomic decomposition by basis pursuit. SIAM Review, 43(1): 129-159.
DOI:10.1137/S003614450037906X |
|
Chopra S, Castagna J P, Xu Y.
2009. Thin-bed reflectivity inversion and some applications. First Break, 27(5): 55-62.
DOI:10.3997/1365-2397.2009009 |
|
Fang Y, Zhang F Q, Wang Y C.
2016. Generalized linear joint PP-PS inversion based on two constraints. Applied Geophysics, 13(1): 103-105, 220.
DOI:10.1007/s11770-016-0527-3 |
|
Huang H D, Wang Y C, Guo F, et al.
2015. Zoeppritz equation-based prestack inversion and its application in fluid identification. Applied Geophysics, 12(2): 199-211.
DOI:10.1007/s11770-015-0483-3 |
|
Li L M, Luo X X, Wang M C, et al.
2010. 3D PP-PS joint inversion method and application in anisotropic medium. Oil Geophysical Prospecting , 45(1): 60-65.
|
|
Liu X J, Yin X Y, Wu G C, et al.
2016. Identification of deep reservoir fluids based on basis pursuit inversion for elastic impedance. Chinese Journal of Geophysics , 59(1): 277-286.
DOI:10.6038/cjg20160123 |
|
Nguyen T, Castagna J P.
2010. High resolution seismic reflectivity inversion. Journal of Seismic Exploration, 19(4): 303-320.
|
|
Partyka G. 2005. Spectral Decomposition, Spring 2005 SEG Distinguished Lecturer. http://ce.seg.org/dl/spring2005/partykaabstract.shtml.
|
|
Portniaguine O, Castagna J. 2005. Spectral inversion:lessons from modeling and Boonesville case study.//75th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 1638-1641, doi:10.1190/1.2148009.
|
|
Puryear C I, Castagna J P.
2008. Layer-thickness determination and stratigraphic interpretation using spectral inversion:Theory and application. Geophysics, 73(2): R37-R48.
DOI:10.1190/1.2838274 |
|
Wei X C, Liu T, Chen T S, et al. 2012. AVO inversion in ray parameter domain.//China Oil and Gas Forum-Geophysical Exploration Technology Symposium. 30-38.
|
|
Widess M B.
1973. How thin is a thin bed?. Geophysics, 38(6): 1176-1180.
DOI:10.1190/1.1440403 |
|
Yin X Y, Liu X J, Wu G C, et al.
2016. Basis pursuit inversion method under model constraint. Geophysical Prospecting for Petroleum, 55(1): 115-122.
DOI:10.3969/j.issn.1000-1441.2016.01.015 |
|
Zhang F Q, Wei F J, Wang Y C, et al.
2013. Generalized linear AVO inversion with the priori constraint of trivariate Cauchy distribution based on Zoeppritz equation. Chinese Journal of Geophysics , 56(6): 2098-2115.
DOI:10.6038/cjg20130630 |
|
Zhang R, Castagna J P.
2011. Seismic sparse-layer reflectivity inversion using basic pursuit decomposition. Geophysics, 76(6): 147-158.
DOI:10.1190/geo2011-0103.1 |
|
Zhang R, Sen M K, Srinivasan S.
2013a. Multi-trace basis pursuit inversion with spatial regularization. Journal of Geophysics and Engineering, 10(3): 35012-35017.
DOI:10.1088/1742-2132/10/3/035012 |
|
Zhang R, Sen M K, Srinivasan S.
2013b. A prestack basis pursuit seismic inversion. Geophysics, 78(1): R1-R11.
DOI:10.1190/geo2011-0502.1 |
|
Zong Z Y, Yin X Y, Zhang F, et al.
2012. Reflection coefficient equation and pre-stack seismic inversion with Young's modulus and Poisson ratio. Chinese Journal of Geophysics , 55(11): 3786-3794.
DOI:10.6038/j.issn.0001-5733.2012.11.025 |
|
李录明, 罗省贤, 王明春.
2010. 各向异性介质三维纵横波联合叠前反演方法及应用. 石油地球物理勘探, 45(1): 60–65.
|
|
刘晓晶, 印兴耀, 吴国忱, 等.
2016. 基于基追踪弹性阻抗反演的深部储层流体识别方法. 地球物理学报, 59(1): 277–286.
DOI:10.6038/cjg20160123 |
|
魏修成, 刘韬, 陈天胜等. 2012. 射线参数域AVO反演方法. //中国油气论坛-地球物理勘探技术专题研讨会. 30-38.
|
|
印兴耀, 刘晓晶, 吴国忱, 等.
2016. 模型约束基追踪反演方法. 石油物探, 55(1): 115–122.
DOI:10.3969/j.issn.1000-1441.2016.01.015 |
|
张丰麒, 魏福吉, 王彦春, 等.
2013. 基于精确Zoeppritz方程三变量柯西分布先验约束的广义线性AVO反演. 地球物理学报, 56(6): 2098–2115.
DOI:10.6038/cjg20130630 |
|
宗兆云, 印兴耀, 张峰, 等.
2012. 杨氏模量和泊松比反射系数近似方程及叠前地震反演. 地球物理学报, 55(11): 3786–3794.
DOI:10.6038/j.issn.0001-5733.2012.11.025 |
|