CN108614295B - 一种基于广义地震子波的地层q值计算方法 - Google Patents

一种基于广义地震子波的地层q值计算方法 Download PDF

Info

Publication number
CN108614295B
CN108614295B CN201810460125.0A CN201810460125A CN108614295B CN 108614295 B CN108614295 B CN 108614295B CN 201810460125 A CN201810460125 A CN 201810460125A CN 108614295 B CN108614295 B CN 108614295B
Authority
CN
China
Prior art keywords
wavelet
frequency
stratum
attenuation
fractional order
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810460125.0A
Other languages
English (en)
Other versions
CN108614295A (zh
Inventor
张金淼
高静怀
翁斌
姜秀娣
桑淑云
张国伟
陈剑军
王清振
丁继才
何洋洋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Offshore Oil Corp CNOOC
CNOOC Research Institute Co Ltd
Original Assignee
China National Offshore Oil Corp CNOOC
CNOOC Research Institute Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China National Offshore Oil Corp CNOOC, CNOOC Research Institute Co Ltd filed Critical China National Offshore Oil Corp CNOOC
Priority to CN201810460125.0A priority Critical patent/CN108614295B/zh
Publication of CN108614295A publication Critical patent/CN108614295A/zh
Application granted granted Critical
Publication of CN108614295B publication Critical patent/CN108614295B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及一种基于广义地震子波的地层Q值计算方法，包括以下步骤：1)对采集到的零偏VSP资料进行预处理，进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；3)计算每一个地层的Q值。本发明基于广义地震子波假设，从粘弹介质中单程波传播理论出发，推导得到了一种利用子波包络峰值处瞬时频率的变化来计算地层Q值的新方法，由于广义地震子波具有参考频率和分数阶系数两个灵活的参数，可以选取最佳的参数来匹配实际地震子波的振幅谱，从而使得估计方法具有良好的子波适应性，并且提高了Q值估计的精度。

Description

一种基于广义地震子波的地层Q值计算方法
技术领域
本发明涉及一种基于广义地震子波的地层Q值计算方法，属于地震勘探技术领域。
背景技术
地震波在地下吸收介质中的传播会经历衰减和频散，这种粘弹性效应一般由品质因子Q来定量表示。在地震勘探中，准确估计Q值是非常重要的。Q可以作为烃类检测和储层描述的工具，也可以作为反Q滤波，非平稳反褶积或逆时偏移中实现吸收补偿的重要参数。
在实际中，Q常常假设为在地震频段内是与频率无关的，也就是所谓的常Q模型。当地震波在吸收介质中传播，其高频分量较低频分量衰减得快。因此，可以通过地震波的频率成分的变化来计算地震衰减，如对数谱比法(Logarithm spectral ratio,LSR)、中心频率偏移法(Centroid frequency shift,CFS)和峰值频率偏移法(Peak frequency shift,PFS)等。
对数谱比法对震源子波没有先验假设，但是其对噪音特别敏感，需要选择合适的带宽范围来拟合斜率。中心频率偏移法假设震源子波振幅谱是高斯形状的，由于其使用幅度谱的统计特征来估计衰减，因而估计结果较稳定，但如果震源子波不能满足高斯形状，那么此方法有较大的系统误差。峰值频率偏移法则假设震源子波是Ricker子波，由于振幅谱峰值位置容易受到噪音影响，因而峰值频率偏移法不够稳健。对CFS和PFS而言，其基于简单的震源子波假设可能导致不适应性问题。
Yang和Gao提出了一种利用包络峰值处瞬时频率(Envelope peak instantaneousfrequency，EPIF)偏移计算Q值的方法，在计算时可以避免加时窗问题，提高了Q值估计的纵向分辨率。然而，此方法的前提是假设子波谱可以用一个常相位的高斯函数表示。实际上，子波谱常常并不是高斯的，而且随着波在吸收介质中的传播，子波谱的形状也会被显著的改变。
综上，需要具有更好的子波适应性的Q值计算方法，从而提高Q值计算精度。
发明内容
针对上述问题，本发明的目的是提供一种基于广义地震子波的地层Q值计算方法，该方法具有更好的子波适应性，从而使得Q值计算精度更高。
为实现上述目的，本发明采用以下技术方案：一种基于广义地震子波的地层Q值计算方法，包括以下步骤：
1)对采集到的零偏VSP资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；
2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；
3)计算每一个地层的Q值，单个地层Q值的计算方法如下：
①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数
Figure BDA0001660695360000021
和参考频率
Figure BDA0001660695360000022
底部直达波的分数阶系数
Figure BDA0001660695360000023
和参考频率
Figure BDA0001660695360000024
②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为
Figure BDA0001660695360000025
和
Figure BDA0001660695360000026
从而地震波在此地层中的传播时间为：
Figure BDA0001660695360000027
③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率或者加权瞬时频率，分别拾取包络峰值处
Figure BDA0001660695360000028
和
Figure BDA0001660695360000029
时刻对应的瞬时频率值，记为
Figure BDA00016606953600000210
和
Figure BDA00016606953600000211
④计算第n个地层的品质因子Qn为：
Figure BDA00016606953600000212
其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：
Figure BDA00016606953600000213
式中的Γ()是Gamma函数；
衰减前后的包络峰值处瞬时频率为
Figure BDA00016606953600000214
均值参考频率
Figure BDA00016606953600000215
和均值分数阶系数
Figure BDA00016606953600000216
分别为顶部和底部相应参数的均值：
Figure BDA00016606953600000217
进一步地，所述步骤3)的④中，Q值的计算公式的建立过程如下：
根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足
Figure BDA00016606953600000218
其中，f是频率；用广义地震子波来表示归一化的震源子波振幅谱：
Figure BDA00016606953600000219
式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率；
广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：
Figure BDA0001660695360000031
其中，Γ()是Gamma函数：
Figure BDA0001660695360000032
式中的x表示函数自变量，t表示积分变量；
衰减后的振幅谱A(f,τ)改写为
Figure BDA0001660695360000033
衰减后子波的包络峰值处瞬时频率等于其幅度谱加权的平均频率
Figure BDA0001660695360000034
在上式中，作变量替换，令
Figure BDA0001660695360000035
则有
Figure BDA0001660695360000036
在上式中，
Figure BDA0001660695360000037
可以认为是一个小量，引入如下的近似
Figure BDA0001660695360000038
Figure BDA0001660695360000039
从而得到
Figure BDA00016606953600000310
Figure BDA0001660695360000041
从而，衰减后子波的包络峰值处瞬时频率为：
Figure BDA0001660695360000042
对上式中的分子分母同时除以
Figure BDA0001660695360000043
并引入近似
Figure BDA0001660695360000044
衰减后子波的包络峰值处瞬时频率可以进一步化简为：
Figure BDA0001660695360000045
注意到衰减前子波的包络峰值处瞬时频率为
Figure BDA0001660695360000046
那么，衰减后子波的包络峰值处瞬时频率可以改写为
Figure BDA0001660695360000047
反过来，可以由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q为：
Figure BDA0001660695360000048
其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：
Figure BDA0001660695360000049
应用时，可用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同，因此，对衰减前后子波的参考频率和分数阶系数求均值：
Figure BDA0001660695360000051
其中，
Figure BDA0001660695360000052
和
Figure BDA0001660695360000053
分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数；从而，地层品质因子Q计算公式修正为：
Figure BDA0001660695360000054
进一步地，所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：
对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2，利用功率谱计算得到子波的均值频率fm和偏差fσ：
Figure BDA0001660695360000055
Figure BDA0001660695360000056
理论上，广义地震子波的均值频率为
Figure BDA0001660695360000057
Figure BDA0001660695360000058
那么
Figure BDA0001660695360000059
其中，与Gamma函数相关的项可以展开为
Figure BDA00016606953600000510
从而
Figure BDA0001660695360000061
这是一个关于
Figure BDA0001660695360000062
的多项式求根问题，可以用迭代法求得最优的u；进一步，可以计算得到广义地震子波的参考频率：
Figure BDA0001660695360000063
对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数
Figure BDA0001660695360000064
和参考频率
Figure BDA0001660695360000065
求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数
Figure BDA0001660695360000066
和参考频率
Figure BDA0001660695360000067
本发明由于采取以上技术方案，其具有以下优点：本发明基于广义地震子波假设，从粘弹介质中单程波传播理论出发，推导得到了一种利用子波包络峰值处瞬时频率的变化来计算地层Q值的新方法，由于广义地震子波具有参考频率和分数阶系数两个灵活的参数，可以选取最佳的参数来匹配实际地震子波的振幅谱，从而使得估计方法具有良好的子波适应性，并且提高了Q值估计的精度。
计算得到的地层Q值，将有利于含气性预测，或者用于衰减效应的补偿。
附图说明
图1是本发明方法的流程示意图；
图2是零偏VSP资料地层划分示意图；
图3是震源子波分数阶系数变化时，分别采用本发明方法、包络峰值频率法、中心频率偏移法和峰值频率偏移法的Q值计算结果的相对误差的变化曲线比较图；
图4是实际零偏VSP资料的下行波与直达波；其中，图(a)是下行波，图(b)是直达波；
图5是本发明方法计算得到的地层Q值曲线。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
如图1所示，本发明提出了一种基于广义地震子波的地层Q值计算方法，包括以下步骤：
1)对采集到的零偏VSP(Vertical Seismic Profiling，垂直地震剖面)资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；
2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层(如图2所示)；
3)计算每一个地层的Q值，单个地层Q值的计算方法如下：
①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数
Figure BDA0001660695360000071
和参考频率
Figure BDA0001660695360000072
底部直达波的分数阶系数
Figure BDA0001660695360000073
和参考频率
Figure BDA0001660695360000074
②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为
Figure BDA0001660695360000075
和
Figure BDA0001660695360000076
从而地震波在此地层中的传播时间为：
Figure BDA0001660695360000077
③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率或者加权瞬时频率，分别拾取包络峰值处
Figure BDA0001660695360000078
和
Figure BDA0001660695360000079
时刻对应的瞬时频率值，记为
Figure BDA00016606953600000710
和
Figure BDA00016606953600000711
④计算第n个地层的品质因子Qn为：
Figure BDA00016606953600000712
其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：
Figure BDA00016606953600000713
式中的Γ()是Gamma函数；
衰减前后的包络峰值处瞬时频率为
Figure BDA00016606953600000714
均值参考频率
Figure BDA00016606953600000715
和均值分数阶系数
Figure BDA00016606953600000716
分别为顶部和底部相应参数的均值：
Figure BDA00016606953600000717
上述实施例中，步骤3)的④中，Q值的计算公式的建立过程如下：
根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足
Figure BDA00016606953600000718
其中，f是频率。假设归一化的震源子波振幅谱可以用广义地震子波来表示：
Figure BDA00016606953600000719
式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率。
根据Barnes，广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：
Figure BDA0001660695360000081
其中，Γ()是Gamma函数：
Figure BDA0001660695360000082
式中的x表示函数自变量，t表示积分变量；
衰减后的振幅谱A(f,τ)改写为
Figure BDA0001660695360000083
衰减后子波的包络峰值处瞬时频率等于其幅度谱加权的平均频率
Figure BDA0001660695360000084
在上式中，作变量替换，令
Figure BDA0001660695360000085
则有
Figure BDA0001660695360000086
在上式中，
Figure BDA0001660695360000087
可以认为是一个小量，引入如下的近似
Figure BDA0001660695360000088
Figure BDA0001660695360000089
从而得到
Figure BDA00016606953600000810
Figure BDA0001660695360000091
从而，衰减后子波的包络峰值处瞬时频率为：
Figure BDA0001660695360000092
对上式中的分子分母同时除以
Figure BDA0001660695360000093
并引入近似
Figure BDA0001660695360000094
衰减后子波的包络峰值处瞬时频率可以进一步化简为：
Figure BDA0001660695360000095
注意到衰减前子波的包络峰值处瞬时频率为
Figure BDA0001660695360000096
那么，衰减后子波的包络峰值处瞬时频率可以改写为
Figure BDA0001660695360000097
反过来，可以由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q为：
Figure BDA0001660695360000098
其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：
Figure BDA0001660695360000099
实际中，可以用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同。因此，可以对衰减前后子波的参考频率和分数阶系数求均值：
Figure BDA0001660695360000101
其中，
Figure BDA0001660695360000102
和
Figure BDA0001660695360000103
分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数。从而，地层品质因子Q计算公式修正为：
Figure BDA0001660695360000104
通过上述公式可知，为了估计地层Q值，需要已知衰减前后的子波包络峰值处的瞬时频率fep(0)和fep(τ)，旅行时间τ，以及用广义地震子波拟合衰减前后子波得到的平均参考频率和平均分数阶系数。
上述实施例中，所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：
对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2。利用功率谱计算得到子波的均值频率fm和偏差fσ：
Figure BDA0001660695360000105
Figure BDA0001660695360000106
理论上，广义地震子波的均值频率应为
Figure BDA0001660695360000107
Figure BDA0001660695360000108
那么
Figure BDA0001660695360000109
其中，与Gamma函数相关的项可以展开为
Figure BDA0001660695360000111
从而
Figure BDA0001660695360000112
这是一个关于
Figure BDA0001660695360000113
的多项式求根问题，可以用迭代法求得最优的u。进一步，可以计算得到广义地震子波的参考频率：
Figure BDA0001660695360000114
利用上述方法，可以分别求得衰减前后的分数阶系数与参考频率。对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数
Figure BDA0001660695360000115
和参考频率
Figure BDA0001660695360000116
求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数
Figure BDA0001660695360000117
和参考频率
Figure BDA0001660695360000118
上述实施例中，步骤3)的②中，计算信号的包络的过程如下：
对于待分析的实地震道x(t)，其对应的复地震道为：
z(t)＝x(t)+iy(t)＝a(t)exp[iθ(t)]
其中，
Figure BDA0001660695360000119
是正交道，
Figure BDA00016606953600001110
表示Hilbert变换，θ(t)是瞬时相位，定义为
Figure BDA00016606953600001111
a(t)是信号的包络，也称为瞬时振幅，定义为
Figure BDA00016606953600001112
根据上述方法，可以分别计算衰减前后子波的包络。对于第n个地层而言，求得的衰减前子波的包络即为地层顶部接收到的信号包络，求得的衰减后子波的包络即为地层底部接收到的信号包络。
上述实施例中，步骤3)的③中，对接收到的信号计算顺时频率、加权顺时频率的过程如下：
对待分析的信号，瞬时频率定义为
Figure BDA0001660695360000121
其中，ε2是阻尼因子，选取为a2(t)最大值的0.00001～0.001倍。为进一步提高瞬时频率的抗噪性能，可以利用包络的平方作为权系数，计算加权瞬时频率：
Figure BDA0001660695360000122
其中W(t)＝a2(t)是加权系数，T是加权窗口的大小。
下面利用合成数据测试本发明方法(GSW)的技术效果：
设地层品质因子为Q＝50，地层旅行时为0.3秒。考虑当震源子波变化时，各种Q值计算方法估计结果的准确性，其中用相对误差来度量准确性。假设震源子波的参考频率为60Hz，分数阶系数从1到2.5变化，变化的步长为0.1，也即有16种情况。分别利用本发明的方法(GSW)，包络峰值频率法(EPIF)，中心频率偏移法(CFS)和峰值频率偏移法(PFS)计算Q值，相对误差(Relative error)随分数阶系数(Fractional value)的变化曲线如图3所示。从图中可以看到，对于EPIF和CFS，它们假设要求震源子波振幅谱是高斯的，因此估计结果存在明显的系统误差。对于PFS，它假设要求振幅谱是Ricker子波，注意到当分数阶系数为u＝2时广义地震子波等价于Ricker子波，此时PFS方法的相对误差很小，但分数阶系数为其他值时，PFS计算结果的误差较大。本发明提出的方法具有最小的相对误差，表明方法的准确性和子波适应性。
本发明的方法用于实际零偏VSP资料的地层Q值计算：
此零偏VSP资料由炸药震源激发，在深度1500-3480米共记录了100道数据，道间距为20米，数据的采样时间步长是1毫秒。图4展示的是实际零偏VSP数据的下行波和拾取的直达波。利用本发明提出的方法计算地层Q值，得到的Q值曲线如图5所示。在3395-3460米范围(图中箭头标记)的低Q值区域，在实际钻井中显示为含气砂岩地层。因此，利用本发明提出的方法计算得到地层Q值，可以帮助解释人员进行含气性预测。
上述各实施例仅用于说明本发明，其中方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (3)

1.一种基于广义地震子波的地层的品质因子Q值计算方法，包括以下步骤：
1)对采集到的零偏VSP资料进行预处理，利用中值滤波进行波场分离得到VSP下行波记录，然后利用窗函数提取每一道的直达波；
2)以井中相邻或者预先设定的检波器的距离，将地下介质划分为多个地层；
3)计算每一个地层的品质因子Q值，单个地层的品质因子Q值的计算方法如下：
①对第n个地层的顶部和底部接收到的直达波，分别利用广义地震子波进行拟合，得到顶部直达波的分数阶系数
Figure FDA0002298399530000011
和参考频率
Figure FDA0002298399530000012
底部直达波的分数阶系数
Figure FDA0002298399530000013
和参考频率
Figure FDA0002298399530000014
②对第n个地层顶部和底部接收到的信号，分别计算其包络，拾取包络峰值处的时间位置，分别为
Figure FDA0002298399530000015
和
Figure FDA0002298399530000016
从而地震波在此地层中的传播时间为：
Figure FDA0002298399530000017
③对第n个地层顶部和底部接收到的信号，分别计算其瞬时频率和幅度谱加权平均频率，分别拾取包络峰值处
Figure FDA0002298399530000018
和
Figure FDA0002298399530000019
时刻对应的瞬时频率值，记为
Figure FDA00022983995300000110
和
Figure FDA00022983995300000111
④计算第n个地层的品质因子Qn为：
Figure FDA00022983995300000112
其中，L(u)是关于广义地震子波分数阶系数的函数，其表达式为：
Figure FDA00022983995300000113
式中的Γ()是Gamma函数；
衰减前后的包络峰值处瞬时频率为
Figure FDA00022983995300000114
均值参考频率
Figure FDA00022983995300000115
和均值分数阶系数
Figure FDA00022983995300000116
分别为顶部和底部相应参数的均值：
Figure FDA00022983995300000117
2.如权利要求1所述的一种基于广义地震子波的地层的品质因子Q值计算方法，其特征在于：所述步骤3)的④中，地层的品质因子Q值的计算公式的建立过程如下：
根据粘弹介质中单程波的传播理论，对于在品质因子为Q的衰减介质中传播的震源子波，设其振幅谱为S(f)，经过时间τ，接收到的子波振幅谱A(f,τ)满足
Figure FDA00022983995300000118
其中，f是频率；用广义地震子波来表示归一化的震源子波振幅谱：
Figure FDA0002298399530000021
式中，u是广义地震子波的分数阶系数；f0是广义地震子波的参考频率；
广义地震子波包络峰值处的瞬时频率等于它的幅度谱加权平均频率：
Figure FDA0002298399530000022
其中，Γ()是Gamma函数：
Figure FDA0002298399530000023
式中的x表示函数自变量，t表示积分变量；
衰减后的振幅谱A(f,τ)改写为
Figure FDA0002298399530000024
衰减后子波的包络峰值处瞬时频率等于其幅度谱加权平均频率
Figure FDA0002298399530000025
在上式中，作变量替换，令
Figure FDA0002298399530000026
则有
Figure FDA0002298399530000027
在上式中，
Figure FDA0002298399530000028
是一个小量，引入如下的近似
Figure FDA0002298399530000029
Figure FDA0002298399530000031
从而得到
Figure FDA0002298399530000032
Figure FDA0002298399530000033
从而，衰减后子波的包络峰值处瞬时频率为：
Figure FDA0002298399530000034
对上式中的分子分母同时除以
Figure FDA0002298399530000035
并引入近似
Figure FDA0002298399530000036
衰减后子波的包络峰值处瞬时频率进一步化简为：
Figure FDA0002298399530000037
注意到衰减前子波的包络峰值处瞬时频率为
Figure FDA0002298399530000038
那么，衰减后子波的包络峰值处瞬时频率改写为
Figure FDA0002298399530000039
反过来，由经过传播时间τ前后，子波包络峰值处的瞬时频率来估计地层的品质因子Q值为：
Figure FDA00022983995300000310
其中，系数函数L(u)是关于广义地震子波分数阶系数的函数，表示为：
Figure FDA0002298399530000041
应用时，用广义地震子波分别来拟合衰减前子波与衰减后的子波，但一般地，衰减前后子波的参考频率和分数阶系数并不相同，因此，对衰减前后子波的参考频率和分数阶系数求均值：
Figure FDA0002298399530000042
其中，
Figure FDA0002298399530000043
和
Figure FDA0002298399530000044
分别是衰减前后子波的参考频率，us和ur分别是衰减前后子波的分数阶系数；从而，地层的品质因子Q值计算公式修正为：
Figure FDA0002298399530000045
3.如权利要求2所述的一种基于广义地震子波的地层的品质因子Q值计算方法，其特征在于：所述步骤3)的①中，利用广义地震子波拟合得到地层顶部直达波和底部直达波的分数阶系数、参考频率的过程如下：
对待拟合的实际地震子波，利用快速傅里叶变换计算其幅度谱为|H(f)|，则功率谱为|H(f)|2，利用功率谱计算得到子波的均值频率fm和偏差fσ：
Figure FDA0002298399530000046
Figure FDA0002298399530000047
理论上，广义地震子波的均值频率为
Figure FDA0002298399530000048
Figure FDA0002298399530000049
那么
Figure FDA00022983995300000410
其中，与Gamma函数相关的项展开为
Figure FDA0002298399530000051
从而
Figure FDA0002298399530000052
这是一个关于
Figure FDA0002298399530000053
的多项式求根问题，用迭代法求得最优的u；进一步计算得到广义地震子波的参考频率：
Figure FDA0002298399530000054
对于第n个地层而言，求得的衰减前的分数阶系数与参考频率即为地层顶部直达波的分数阶系数
Figure FDA0002298399530000055
和参考频率
Figure FDA0002298399530000056
求得的衰减后的分数阶系数与参考频率即为地层底部直达波的分数阶系数
Figure FDA0002298399530000057
和参考频率
Figure FDA0002298399530000058
CN201810460125.0A 2018-05-15 2018-05-15 一种基于广义地震子波的地层q值计算方法 Active CN108614295B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810460125.0A CN108614295B (zh) 2018-05-15 2018-05-15 一种基于广义地震子波的地层q值计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810460125.0A CN108614295B (zh) 2018-05-15 2018-05-15 一种基于广义地震子波的地层q值计算方法

Publications (2)

Publication Number Publication Date
CN108614295A CN108614295A (zh) 2018-10-02
CN108614295B true CN108614295B (zh) 2020-04-21

Family

ID=63663228

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810460125.0A Active CN108614295B (zh) 2018-05-15 2018-05-15 一种基于广义地震子波的地层q值计算方法

Country Status (1)

Country Link
CN (1) CN108614295B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109188524A (zh) * 2018-10-26 2019-01-11 辽宁工程技术大学 考虑介质品质因子的地震中场地反应计算方法
CN109188523B (zh) * 2018-10-26 2020-03-17 辽宁工程技术大学 考虑介质对地震波吸收衰减的地震中场地反应计算方法
CN111427080B (zh) * 2020-03-13 2021-08-13 王仰华 地震数据空变广义子波的提取方法
CN111427081B (zh) * 2020-03-28 2022-04-08 王仰华 时变地震衰减模型的构建方法
CN114252910A (zh) * 2020-09-25 2022-03-29 中国石油天然气股份有限公司 地震品质因子确定方法及装置
CN113703049B (zh) * 2021-08-27 2024-04-23 中铁二十局集团安哥拉国际有限责任公司 一种基于n次傅里叶谱的子波估计方法
CN114152985B (zh) * 2021-12-14 2022-10-28 中国地质大学（北京） 一种确定地下古河道边界及其内部薄砂体厚度的方法
CN114563824B (zh) * 2022-02-25 2024-01-30 成都理工大学 一种二阶多重同步挤压多项式chirplet变换薄储层识别方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102323615B (zh) * 2011-06-02 2013-11-06 中国石油集团川庆钻探工程有限公司地球物理勘探公司 利用地震数据进行储层预测和流体识别的方法和装置
CN102221708B (zh) * 2011-06-03 2013-05-08 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于分数阶傅里叶变换的随机噪声压制方法
US9268047B2 (en) * 2012-10-12 2016-02-23 Rock Solid Images, Inc Geophysical surveying
CN102928875B (zh) * 2012-11-07 2016-05-11 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于分数阶傅里叶域的子波提取方法
CN103728662A (zh) * 2014-01-03 2014-04-16 中国海洋石油总公司 一种基于地震信号包络峰值的地层介质品质因子估计方法
US20150355358A1 (en) * 2014-06-06 2015-12-10 Schlumberger Technology Corporation Generalized spectral decomposition
CN105388523A (zh) * 2014-09-04 2016-03-09 中国石油化工股份有限公司 一种高精度的品质因子提取方法

Also Published As

Publication number Publication date
CN108614295A (zh) 2018-10-02

Similar Documents

Publication Publication Date Title
CN108614295B (zh) 一种基于广义地震子波的地层q值计算方法
CN103376464B (zh) 一种地层品质因子反演方法
US8553498B2 (en) Q tomography method
CN109254324B (zh) 全频保幅地震数据处理方法和装置
CN109669212B (zh) 地震数据处理方法、地层品质因子估算方法与装置
CN104237945B (zh) 一种地震资料自适应高分辨处理方法
CN101923176B (zh) 一种利用地震数据瞬时频率属性进行油气检测的方法
CN105044777B (zh) 基于经验模态分解检测地震标志层强反射振幅消除的方法
CN110967749A (zh) 一种vsp地震资料频变q值估计与反q滤波方法
CN103592680B (zh) 一种基于正反演的测井数据和深度域地震剖面合成方法
Wang et al. Improving the resolution of seismic traces based on the secondary time–frequency spectrum
CN104808245A (zh) 道集优化处理方法及其装置
CN105044772A (zh) 一种实现时变谱模拟反褶积的方法和装置
CN109031414B (zh) 一种基于希尔伯特变换的振幅增益方法及处理终端
Zhang et al. Interval Q inversion based on zero-offset VSP data and applications
Chen et al. Improving the Precision of Surface Seismic Data Processing by Walkaway VSP
Wang et al. A method for absorption compensation based on adaptive molecular decomposition
CN112684501B (zh) 一种基于谱比面积的q值估计方法与应用
CN110749923A (zh) 一种基于范数方程提高分辨率的反褶积方法
CN110568491B (zh) 一种品质因子q的估算方法
CN114137606A (zh) 一种稳健的谱模拟反褶积方法
CN112241025B (zh) 一种井震联合地层压力确定方法及系统
CN112526611A (zh) 表层地震波品质因子的提取方法及装置
Pramanik et al. Estimation of Q from borehole data and its application to enhance surface seismic resolution: A case study
CN117724165B (zh) 一种基于时变子波的品质因子估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant