2017, Vol. 60
2. 中国石油大学 (北京) CNPC物探重点实验室, 北京 102249;
3. 中国石油大学 (北京) 地球物理与信息工程学院, 北京 102249
2. CNPC Key Lab of Geophysical Exploration, China University of Petroleum, Beijing 102249, China;
3. College of Geophysics and Information Engineering, China University of Petroleum, Beijing 102249, China
近年来，随着石油勘探程度的不断提高，构造油气藏已逐渐开发殆尽，地球物理勘探被迫从寻找构造油气藏转向寻找岩性油气藏，地震反演也逐渐从叠后转向叠前.叠前资料中不仅包含着反射振幅随角度变化的趋势，还携带了更丰富的地下地质信息，包括纵波速度、横波速度和密度等，从而有利于岩性预测和流体识别.
相比于叠后，虽然叠前数据中所包含的反射振幅随角度变化的趋势为反演提供了更丰富的信息，但受到子波带宽的限制，反演仍面临着信息量不足的问题，导致了结果的非唯一性.解决该问题的常用方式是加入先验信息，以约束反演取得期望的结果.例如，Downton和Lines (2004)、Theune等 (2010)、Alemie和Sacchi (2011)、杨培杰等 (2009)假设叠前参数符合柯西分布或改进的柯西分布，从而获得稀疏的结果，提高了反演分辨率.Buland和More (2003)假设叠前参数的先验分布及似然函数均为高斯分布，由此得到参数的后验概率也为高斯分布，并推导了均值和协方差矩阵，进而给出了参数在95%置信水平下的置信区间，从概率论的角度对结果进行了分析.与多数叠前反演方法不同，Ma (2002)直接反演纵横波阻抗而非其变化率，为提高反演的稳定性，Ma在目标函数中加入了纵横波阻抗的低频背景信息作为约束.Sen和Stoffa (1991)通过对地下介质进行等厚度划分的方式进行反演，即需要指定界面个数.Pérez等 (2013)虽没有上述严格的限制，但仍假设反射系数是稀疏的，从而求得叠前参数的位置和幅值.对于反演而言，尽可能地利用来自于测井资料或岩石物理实验的先验信息，通过数学的方式加入一些先验约束固然可以缩小解的范围，指导反演获得期望的结果，但在实际地震资料处理中，这些先验信息往往难以提供或满足，放松假设条件单纯地利用地震数据进行反演不仅适用于更多的情形，还可以作为质量控制的工具.
除了引入先验信息以约束解的范围，采用的算法也是反演取得成功的关键.由于反射界面的位置未知，叠前反演本质上属于非线性反演.Alemie和Sacchi (2011)通过迭代更新的方式进行求解，具有较快的收敛速度.印兴耀等 (1994)采用人工神经网络进行反演，Sen和Stoffa (1991)以及Ma (2002)采用模拟退火算法求得全局最优解，由于采用了这些非线性优化算法，目标函数中更容易加入各种先验约束.印兴耀等 (2016)基于逆算子估计的方法进行AVO反演，相比于一般的优化类算法，具有较高的运算效率.
匹配滤波方法是数字信号处理领域广泛采用的处理技术，其准则是使得输出信号的信噪比达到最大，常用于判断输入信号中是否存在某一特殊波形.Li和Stoica (1996)进一步提出逐点匹配滤波的思想，在对信号的谱进行估计的过程中，针对每个频率设计一个滤波器，从而将非线性问题化解为多个线性问题.Guo和Wang (2012)将该方法用于地震叠后反演，假设每一采样时刻都可能存在反射界面并计算其幅值，理论推导和实验结果均表明该方法在真实界面位置可以得到近似无偏的结果.本文基于逐点匹配滤波思想，对每个采样时刻进行叠前反演并计算残差值，通过在残差曲线上自动拾取波谷的方式确定反射界面位置，然后在频率域采用最小二乘反演对叠前参数进行更准确地估计，从而将非线性问题转化为两步线性反演实现，不仅具有较高的计算效率、较强的稳定性，并且不需要界面个数、待反演参数的概率分布等先验约束，正则化参数的选择对于结果的影响也相对较小，从而有利于反演的实现.
本文首先对逐点AVA反演方法进行了理论推导，明确了残差曲线上波谷位置与反射界面之间的关系，给出了反射界面自动拾取的办法，推导了频率域最小二乘反演的公式，并详细介绍了“界面定位-参数估计”的两步反演算法.在模型测试部分不仅考虑了抗噪性，还分析了界面个数拾取不准确对反演结果的影响，最终利用实际地震资料验证了该方法的有效性.
2 理论 2.1 频率域逐点AVA反演基于层状介质假设，地震记录s (t)可以表示为反射系数序列r (t)与子波w (t)的褶积，公式为
|
(1) |
在频率域为
|
(2) |
其中，S (ω)、W (ω)、R (ω)分别表示地震记录的谱、子波的谱以及反射系数序列的谱.
假设子波或子波的谱已知，通过谱除的方式可以得到反射系数序列在有限频带内的谱.结合傅里叶变换可得：
|
(3) |
其中，
基于水平层状介质假设，将公式 (3) 扩展至叠前以包含角度信息，并考虑加性噪声的影响，可得公式为
|
(4) |
其中，M表示道数；rk(θm)表示第k个界面对应的角度为θm的反射系数.虽然rk(θ1), …, rk(θM) 均为未知量，但其对应的界面位置tk相同，因而相比于叠后谱反演公式 (3)，公式 (4) 增加了信息量.
基于逐点匹配滤波的思想，在对第k个界面进行处理时，仅仅将该界面的反射视作有效信号，于是有：
|
(5) |
其中，
将公式 (5) 重组，公式为
|
(6) |
其中，
(6) 式可简写为
|
(7) |
其中，d=(R1R2…RM)T，Ak为Ek组成的稀疏块状矩阵，rk=[rk1, …, rkM]T.
根据Shuey方程 (Shuey，1985)，有公式为
|
(8) |
其中，pk、gk分别表示第k个界面对应的截距与梯度.基于公式 (8)，有公式为
|
(9) |
其中，xk=[pk, gk]T，
|
将公式 (9) 代入 (7) 中，得到：
|
(10) |
其中，Gk=AkB.
公式 (10) 建立了模型与数据之间的关系，可以采用匹配滤波方法进行求解.匹配滤波的核心思想在于压制干扰突出有效信号，以使得输出信号的信噪比达到最大，若将干扰的协方差矩阵记为Qk，则第k个界面对应的滤波器为 (Li and Stoica, 1996)：
|
(11) |
其中，H表示复共轭转置.进而求得截距与梯度参数，公式为
|
(12) |
干扰的协方差矩阵通常是利用观测数据自适应地进行估计，根据其估算方式的不同，滤波器分为Capon滤波器、APES滤波器、Fourier滤波器等 (Stoica and Moses, 1997)，本文采用相对简单的滤波器--Fourier滤波器，即将干扰的协方差矩阵近似为单位斜对角阵，从而有：
|
(13) |
因为Gk=AkB，所以GkHGk=BTAkHAkB，而AkHAk=NI(I为单位斜对角阵)，于是GkHGk=NBTB，因而有公式为
|
(14) |
为稳定反演，在公式 (14) 中引入正则化项，于是有公式为
|
(15) |
由于反射界面个数与位置均未知，因而假设每个采样点处都可能存在反射界面，根据采样时刻计算Ak，以利用公式 (15) 估计该采样点对应的截距与梯度向量
公式 (15) 中，有：
|
(16) |
公式 (16) 表明，AkHd/N中的第m(m=1, …, M) 个元素，就是角道集中第m道的谱的反傅里叶变换，因而根据公式 (15) 计算的
进一步将公式 (3) 代入公式 (16) 中，得公式为
|
(17) |
显然，求解向量AkHd/N中的第m个元素
鉴于逐点AVA反演方法的如上不足，进一步利用合成数据与观测数据的残差对结果进行校正.由于该方法针对每个采样时刻进行参数估计，因而每个采样点都对应一个残差值，该值是对反演结果进行评估的重要依据.第k个采样点对应的残差值为
|
(18) |
将公式 (15) 代入公式 (18) 中，进一步整理，得公式为
|
(19) |
其中，
基于此前的分析，若在第k个采样点的位置存在一个真实界面，则AkHd中的各个分量将取得局部极大值，因而 (AkHd)HW(AkHd) 取得局部极大值，由此ek将取得局部极小值.也就是说，基于公式 (19) 对所有采样点计算ek，得到残差e随深度的变化曲线，该曲线上局部极小值即波谷的位置，与反射界面存在着对应关系.并且，反射越强，(AkHd)HW(AkHd) 的取值越大，导致ek越小，因而对逐点反演方法获得的残差曲线进行分析，不仅可以确定界面位置还有助于对界面的反射强度进行判断.虽然基于公式 (20) 计算的残差曲线在真实界面位置附近也会呈现出一定的旁瓣，但由于残差值均大于零，故残差曲线上不存在符号相反的情形，因而界面位置更容易识别.
事实上，真实反射界面位置之所以与残差曲线上的波谷位置相对应，还可以从另一个角度进行解释.对于反演而言，残差大小与观测数据的信噪比紧密相关.信噪比较高，意味着噪声较弱，合成数据与观测数据能更好地进行匹配，从而得到较小的残差值；相反，若信噪比较低，则意味着噪声较强，导致较大的残差值.逐点AVA反演方法的基本思想是，在对某一采样时刻进行参数估计时，仅将该位置处的反射视作有效信号，而将所有其余界面的反射以及噪声的影响均视作干扰.根据公式 (10)，对各个采样点进行反演所基于的数据d是相同的, 因而对于那些反射较弱甚至不存在反射的采样点而言，意味着干扰较强，甚至数据d完全被视作干扰，导致反演的残差较大；而对于那些存在反射尤其是强反射的采样点而言，意味着干扰相对较弱，从而得到较小的残差值.总体而言，反射越弱，相应位置的残差值越大；反射越强，相应位置的残差值越小.反言之，残差值较小，相应位置可能存在较强的反射，而残差值较大甚至接近于‖d‖2，则表明相应位置的反射较弱甚至不存在反射.因而利用残差值大小可以对一个采样点处是否存在反射以及反射的强弱进行分析.
2.3 基于反射界面自动拾取的叠前反演算法由于在逐点AVA反演获得的残差曲线上，波谷与真实界面位置存在着对应关系，因而通过自动拾取波谷的方式可以定位反射界面位置，进而利用最小二乘反演对截距和梯度参数进行更准确的估计.在界面位置已知的情形下，对公式 (4) 进行重组，得到频率域最小二乘反演的公式为
|
(20) |
简写为
|
(21) |
其中，P为给定的界面个数，F=CD，C=
|
(22) |
于是，反演基于两步进行：首先在逐点反演获得的残差曲线上自动拾取波谷以定位反射界面；然后在界面位置已知的情形下进行频率域最小二乘反演.具体步骤为：
(1) 针对所有采样时刻，基于公式 (19) 计算残差，得到残差曲线e.
(2) 在残差曲线上拾取波谷，并按照残差值从小到大的顺序对波谷位置进行排序，得到t1, t2, …, tL，L为拾取的波谷数目.
(3) 选取t1、t2、…、tP(P=1, 2, …, L) 作为界面位置，基于公式 (22) 进行反演，并计算最小二乘反演的残差
(4) 选取残差趋于平稳的临界点，即为最优的界面个数，将其对应的反演结果作为最终结果.
从上述步骤 (1) 可以看出，虽然反射界面的自动拾取是基于逐点AVA反演方法提出的，但却不必先进行AVA反演，只需利用公式 (19) 直接计算残差值，从而节省了计算时间，而且在计算残差的过程中，各个采样点之间相互独立，不存在依赖关系，因而步骤 (1) 可并行完成.
虽然反射界面位置与逐点反演获得的残差曲线上的波谷位置相对应，但受到薄层调谐效应的影响，在残差曲线上两个相距较近的界面之间仍有可能出现“虚假”的波谷，也就是说在这些波谷位置，并不存在真实的反射界面.因而上述步骤 (2) 至 (4) 采用尝试求解，以确定界面个数.根据2.2节的分析，残差值较小，则相应位置可能存在较强的反射，因而步骤 (2) 中基于残差值由小到大的顺序对波谷位置进行排序，并通过残差曲线收敛的方式判断界面个数.由于残差曲线上拾取的波谷数目远远小于采样点数，因而尝试求解的计算成本很低.
此外，需要注意的是，无论是对于逐点反演还是最小二乘反演，正则化参数的选择都非常重要，不同的正则化参数会导致不同的反演结果.该参数通常基于观测数据的信噪比进行选取，若信噪比较低，则选取较大的正则化参数，使得反演更加依赖于先验约束；若信噪比较高，则选取较小的正则化参数，使得反演更加信任观测数据.具体的选择方法有L曲线、Χ2分布等多种方法 (Hansen，1992；Hennenfent et al., 2008；Van Den Berg and Friedlander, 2008)，本文中最小二乘反演的正则化参数λ2的选取采用了L曲线法.对于逐点反演方法而言，虽然不同的正则化参数λ1会导致反演结果不同，从而导致各采样点对应的残差值发生变化，但只要对于所有的采样点给定相同的λ1，就不会影响残差曲线上波谷位置的拾取.因为，根据公式 (19)，相同的λ1意味着所有采样点对应的矩阵W是相同的，增加λ1会导致所有采样点对应的残差值增加，缩小λ1也会导致所有采样点对应的残差值减小，从而不影响波谷位置的拾取.而且，根据2.2节对于AkHd的分析，如果第k个采样点的位置存在一个真实反射界面，则AkHd取得局部极大值，因而根据公式 (19)，残差e随正则化参数λ1单调增加或减小的趋势在真实界面位置表现得更为明显.这一现象也可以从另外一个角度进行解释，由于逐点反演所基于的数据d是相同的，因而对于那些存在反射尤其是强反射的采样时刻，干扰相对较弱，调整正则化参数可以在很大程度上影响反演结果，从而改变合成数据与观测数据的匹配程度，导致残差值受正则化参数的影响较大；相反，对于那些反射较弱甚至不存在反射的采样点而言，干扰较强，甚至整个数据d都成为干扰，在此情况下即使调整正则化参数，也难以改变合成数据与观测数据的匹配程度，导致残差值随正则化参数的变化较小.正因为残差值随正则化参数λ1的变化在真实界面位置表现得更为明显，因而利用这一现象也有助于对界面位置进行识别.
3 数值例子
参照Pérez等 (2013)的模型对本文提出的方法进行了实验分析.图 1a和1b给出了截距和梯度模型，基于主频为30 Hz的零相位雷克子波，利用Shuey近似合成的无噪角道集如图 1c所示，加入标准差为
|
图 1 真实模型与合成数据 (a) 截距；(b) 梯度；(c) 无噪叠前道集；(d)-(f) 不同信噪比的叠前道集. Fig. 1 True model and synthetic data (a) Intercept; (b) Gradient; (c) Noise-free prestack gather; (d)-(f) Prestack gathers with different S/N. |
图 2以S/N=20的角道集为例，展示了反演的过程.对于本文提出的方法，要拾取波谷位置，只需根据公式 (19) 逐点计算残差即可，不必先利用公式 (15) 反演截距和梯度，但为了进行对比，图 2a仍给出了逐点反演的结果 (给定λ1=1)；图 2b中蓝线和红线分别给出了λ1=1和λ1=20的情况下逐点计算的残差曲线，以展示正则化参数的影响；图 2c以λ1=1的残差曲线为例展示了波谷位置的拾取；图 2e和2f分别给出了基于不同界面个数进行最小二乘反演获得的残差曲线以及反演的截距和梯度.对于图 2a所示的逐点反演结果，虽然其波峰和波谷与真实反射系数之间存在一定的对应关系，但正如2.2节理论分析所预测，该结果旁瓣较大，分辨率不足，精度较低.而在图 2b所示的残差曲线上，波谷与真实反射系数之间的对应关系更加明显，并且正则化参数越大，合成数据与观测数据的残差越大 (红线所表示的残差值大于蓝线所表示的残差值)，但二者之间的差异只有在真实反射界面处才较为明显，在其余位置，不同λ1对应的残差均趋向于数据能量‖d‖2而几乎重合，这一特征也有利于反射界面的识别.但是，受薄层调谐效应的影响，残差曲线在真实界面附近仍呈现出一定的旁瓣，虽然较强的 (残差值较小) 波谷与反射界面存在很好的对应关系，但在其余不存在反射的位置，也出现了较弱的 (残差值较大) 波谷，例如图 2c中第二和第三个界面之间的波谷甚至比第一个界面处的波谷更为明显.按照残差值由小到大的顺序对拾取的波谷进行排序，并选取前面P个波谷位置进行最小二乘反演，在图 2d所示的残差f随界面个数P变化的曲线上，显然P=14是曲线收敛的临界点，但为了对结果进行更详细地分析，图 2e展示了选定界面数分别为14、16、18和20反演的截距与梯度.一方面，随着界面个数误差的不断增加，反演的准确度逐渐降低，尤其是梯度的反演结果，会呈现出较多的虚假反射；但另一方面，即使界面个数误差较大，反演的截距却较为稳定和准确，从而有利于反演结果的解释.事实上，关于梯度反演的不稳定性，已有多位学者进行过阐述 (Herrmann and Cambois, 2001；Whitcombe et al., 2004)，也提出了一些解决办法 (Downton and Lines, 2001，2004；Buland and More, 2003)，但该问题并非本文研究的重点.
|
图 2 反演的中间结果 (a) 逐点反演的截距与梯度；(b) λ1=1(蓝线) 和λ1=20(红线) 对应的逐点反演残差曲线; (c) 在λ1=1对应的残差曲线上拾取的波谷位置 (红圈)；(d) 最小二乘反演的残差f随界面个数的变化；(e) 给定不同界面数反演的截距 (第一行) 和梯度 (第二行).(b) 与 (c) 中绿色直线表示数据能量‖d‖2，黑色虚线指明了真实反射系数的位置. Fig. 2 Intermediate results of inversion (a) Intercepts and gradients of sample-by-sample inversion; (b) Residual curves of sample-by-sample inversion with λ1=1 (blue curve) and λ1=20 (red curve); (c) Locations of troughs picked up on the residual curve with λ1=1; (d) Variation of residual f of least-squares inversion with the number of interfaces; (e) Intercepts and gradients (first row and second row, respectively) inverted in the case of different numbers of interfaces. In (b) and (c), green lines denote the energy ‖d‖2 of the observed data and black dash lines show the locations of true interfaces. |
图 3针对S/N=5和S/N=10的叠前道集进行了反演分析.图 3a显示了逐点反演获得的残差曲线，由于信噪比为5的道集中噪声较强，导致了其数据残差较大，即便如此，自动拾取的波谷位置仍与真实反射界面匹配良好.图 3b显示了最小二乘反演的残差曲线，与图 3a类似，信噪比为5的道集对应的残差较大.图 3c给出了不同信噪比的地震数据对应的最终反演结果.信噪比为20和10的叠前道集都得到了近乎完美的结果，对于信噪比为5的道集而言，由于噪声过强，导致梯度反演结果不够精确，但截距反演结果较为准确，在信噪比仅为5的情形下，该结果是可以接受的.
|
图 3 不同信噪比的叠前道集对应的反演结果对比 (a) λ1=1的情形下逐点反演的残差曲线；(b) 最小二乘反演的残差曲线；(c) 反演的截距 (第一行) 与梯度 (第二行). Fig. 3 Comparison of inversion results for prestack gathers with different S/N (a) Residual curves and troughs of sample-by-sample inversion in the case of λ1=1; (b) Residual curves for least-squares inversion; (c) Inverted intercepts and gradients inverted (first row and second row, respectively). |
为验证本方法的实际应用效果，我们对中国大陆陆地某勘探工区的叠前数据进行了处理.图 4显示了中心角度分别为5°、15°和25°的部分叠加剖面，均为131道，采样间隔2 ms.采用自相关法分别从三个剖面提取子波，反演采用的频带为5~65 Hz，得到如图 5所示的反演结果.虽然梯度剖面的连续性及稳定性相对稍差，但反演的截距剖面却展示了较好的横向连续性和较高的分辨率.此外，根据测井资料提供的纵横波速度和密度，计算了第71道所处位置的截距曲线和梯度曲线作为比对标准.首先，反演得到的界面位置与测井曲线所展现的界面位置对应良好，表明了基于残差曲线的波谷进行反射界面自动拾取方法的可行性.其次，反演的截距的极性和幅值与测井资料的计算结果匹配良好，也进一步验证了界面位置拾取的准确性，并表明了频率域最小二乘反演方法的可靠性.
|
图 4 不同角度的地震剖面 (a) 小角度叠加剖面；(b) 中角度叠加剖面；(c) 大角度叠加剖面. Fig. 4 Seismic profiles corresponding to three angles (a) Stacked profile with small angle; (b) Stacked profile with middle angle; (c) Stacked profile with large angle. |
|
图 5 反演结果 (a) 截距剖面；(b) 梯度剖面.黑色曲线为根据测井资料计算的截距和梯度. Fig. 5 Inversion results (a) Intercept profile; (b) Gradient profile. Black curves denote intercept and the gradient calculated from well logs, respectively. |
本文将逐点匹配滤波思想应用于叠前AVA反演，提出了在残差曲线上拾取波谷进行反射界面自动定位的方法，从而将AVA非线性反演转化为两步线性反演.相比于传统的同时估计反射位置和幅值的非线性方法，该方法不需要关于参数的个数、分布等先验假设，并且具有运算速度快、求解稳定等优点.首先，由于针对每个采样点进行的是线性反演，因而计算速度快，并且各个采样点之间可以采用并行处理以进一步提高计算速度.其次，针对每个采样点给定相同的正则化参数，虽然会造成所有采样点对应的残差值增加或减小，但不会影响残差曲线上波谷位置的拾取，从而降低了反演对于正则化参数的依赖，即使是在信噪比较低的情形下，仍获得了较好的反演结果.再有，虽然薄层的调谐效应可能导致拾取的界面个数多于实际的界面个数，但截距的反演结果受界面个数不准的影响较小，从而有利于反演结果的分析与解释.
尽管模型反演与实际资料处理都展示了新方法的有效性，但该方法仍可以进一步完善.一方面，在逐点反演时采用Capon、APES等滤波器可以提高单点反演的准确度，进而提高分辨率；另一方面，对过于复杂的地质模型，可以将逐点反演获得的界面位置作为初值，采用非线性反演提高结果的准确程度.
| Alemie W, Sacchi M D. 2011. High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution. Geophysics, 76(3): R43-R55. DOI:10.1190/1.3554627 | |
| Buland A, Omre H. 2003. Bayesian linearized AVO inversion. Geophysics, 68(1): 185-198. DOI:10.1190/1.1543206 | |
| Downton J E, Lines L R. 2001. Constrained three parameter AVO inversion and uncertainty analysis.//71st SEG Annual Meeting. San Antonio, Texas:SEG, 251-254. | |
| Downton J E, Lines L R. 2004. Three term AVO waveform inversion.//74th SEG Annual Meeting. Denver, Colorado:SEG, 215-218. | |
| Guo R, Wang S X. 2012. A spectral method for reflectivity estimation. J. Geophys. Eng., 9(6): 681-690. DOI:10.1088/1742-2132/9/6/681 | |
| Hansen P C. 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev., 34(4): 561-580. DOI:10.1137/1034115 | |
| Hennenfent G, van den Berg E, Friedlander M P, et al. 2008. New insights into one-norm solvers from the Pareto curve. Geophysics, 73(4): A23-A26. DOI:10.1190/1.2944169 | |
| Herrmann P, Cambois G. 2001. Statistical properties of seismically derived AVO attributes.//71st SEG Annual Meeting. San Antonio, Texas:SEG, 231-234. | |
| Li J, Stoica P. 1996. An adaptive filtering approach to spectral estimation and SAR imaging. IEEE Trans. Signal Process., 44(6): 1469-1484. DOI:10.1109/78.506612 | |
| Ma X Q. 2002. Simultaneous inversion of prestack seismic data for rock properties using simulated annealing. Geophysics, 67(6): 1877-1885. DOI:10.1190/1.1527087 | |
| Pérez D O, Velis D R, Sacchi M D. 2013. High-resolution prestack seismic inversion using a hybrid FISTA least-squares strategy. Geophysics, 78(5): R185-R195. DOI:10.1190/geo2013-0077.1 | |
| Sen M K, Stoffa P L. 1991. Nonlinear one-dimensional seismic waveform inversion using simulated annealing. Geophysics, 56(10): 1624-1638. DOI:10.1190/1.1442973 | |
| Shuey R T. 1985. A simplification of the Zoeppritz equations. Geophysics, 50(4): 609-614. DOI:10.1190/1.1441936 | |
| Stoica P, Moses R L. 1997. Introduction to Spectral Analysis. Upper Saddle River, NJ:Prentice-Hall. | |
| Theune U, Jensås I Ø, Eidsvik J. 2010. Analysis of prior models for a blocky inversion of seismic AVA data. Geophysics, 75(3): C25-C35. DOI:10.1190/1.3427538 | |
| Van Den Berg E, Friedlander M P. 2008. Probing the Pareto frontier for basis pursuit solutions. SIAM J. Sci. Comput., 31(2): 890-912. | |
| Whitcombe D N, Dyce M, McKenzie C J S, et al. 2004. Stabilizing the AVO gradient.//74th SEG Annual Meeting. Denver, Colorado:SEG, 232-235. | |
| Yang P J, Mu X, Yin X Y. 2009. Prestack three-term simultaneous inversion method and its application. Acta Petrolei Sinica (in Chinese), 30(2): 232-236. | |
| Yin X Y, Deng W, Zong Z Y. 2016. AVO inversion based on inverse operator estimation. Chinese J. Geophys. (in Chinese), 59(4): 1457-1468. DOI:10.6038/cjg20160426 | |
| Yin X Y, Wu G C, Zhang H Z. 1994. The application of neural networks in the reservoir prediction. Journal of the University of Petroleum, China (in Chinese), 18(5): 20-26. | |
| 杨培杰, 穆星, 印兴耀. 2009. 叠前三参数同步反演方法及其应用. 石油学报, 30(2): 232–236. | |
| 印兴耀, 邓炜, 宗兆云. 2016. 基于逆算子估计的AVO反演方法研究. 地球物理学报, 59(4): 1457–1468. DOI:10.6038/cjg20160426 | |
| 印兴耀, 吴国忱, 张洪宙. 1994. 神经网络在储层横向预测中的应用. 石油大学学报:自然科学版, 18(5): 20–26. | |