CN110837119B - 一种增强反q值补偿稳定性的地震资料处理方法及系统 - Google Patents

一种增强反q值补偿稳定性的地震资料处理方法及系统 Download PDF

Info

Publication number
CN110837119B
CN110837119B CN201810942502.4A CN201810942502A CN110837119B CN 110837119 B CN110837119 B CN 110837119B CN 201810942502 A CN201810942502 A CN 201810942502A CN 110837119 B CN110837119 B CN 110837119B
Authority
CN
China
Prior art keywords
seismic data
value
reflection coefficient
norm
time
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
CN201810942502.4A
Other languages
English (en)
Other versions
CN110837119A (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 Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
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 Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201810942502.4A priority Critical patent/CN110837119B/zh
Publication of CN110837119A publication Critical patent/CN110837119A/zh
Application granted granted Critical
Publication of CN110837119B publication Critical patent/CN110837119B/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. analysis, for interpretation, for correction
    • 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
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • 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. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G01V1/368—Inverse filtering

Abstract

本发明公开了一种增强反Q值补偿稳定性的地震资料处理方法及系统，包括：步骤1：获取地震数据、反射系数和Q值之间的关系表达式；步骤2：获得Lp范数约束下的反演目标函数；步骤3：获得所述Lp范数约束下的反演目标函数的Q值；步骤4：基于所述关系表达式和所述Lp范数约束下的反演目标函数的Q值，获得反射系数。本发明通过获取地震数据、反射系数和Q值之间的关系表达式以及Lp范数约束下的反演目标函数，计算获取Lp范数约束下的反演目标函数的Q值，基于地震数据、反射系数和Q值之间的关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数，增加Lp范数约束的反演方法，较大地提高了地震信号的分辨率和稳定性好。

Description

一种增强反Q值补偿稳定性的地震资料处理方法及系统
技术领域
本发明属于石油地球物理勘探领域，具体涉及一种增强反Q值补偿稳定性的地震资料处理方法及系统。
背景技术
地震波在地下传播过程中由于地层介质的吸收衰减作用造成能量损失，且高频的衰减比低频严重，导致地震波主频降低，频带变窄，严重影响了中深层资料的成像精度，且在组成物质松散的地层或裂隙发育的地层中，地震波的吸收响应要比地震波速响应更为敏感。介质粘弹性引起的地震波衰减特性通常用品质因子Q值来描述，它与介质内部的结构特征以及饱和度、孔隙度、渗透率等因素密切相关。
地层的Q值滤波效应同地震子波的滤波效应均能降低地震资料的分辨率，导致高频成分的丢失，使得深部反射波能量减弱。因此从未作振幅补偿和相位校正的地震剖面上识别深部同相轴及进一步反演工作均存在困难。反Q值滤波方法是一种消除地层Q值滤波效应的普遍方法，但由于其在振幅补偿方面本质上不稳定，现阶段依然有许多科研人员在解决该问题上攻坚克难。
地震波在传播过程中其振幅呈指数衰减，在传播到深部位置时，有效信号的能量便淹没于背景噪声之中，而一般的反Q滤波方法是对整体能量的指数抬升，因此深部噪声的能量得到急剧放大，从而导致不稳定现象。即便在无噪环境下，由于计算精度问题，数值误差的产生及指数型放大也会带入不稳定现象。近些年来，地球物理学家们提出过各种各样的Q值计算方法，比如谱比法、脉冲振幅法、斜率法和小波域谱比法等，这些方法多数是考虑地震波内部的衰减，它们依然存在计算准确性和稳定性问题。对此特别需要一种合理的反Q值计算方法，能提高反Q值计算的分辨率及稳定性。
发明内容
本发明的目的是提出一种提高地震信号分辨率及稳定性的增强反Q值补偿稳定性的地震资料处理方法及系统。
根据本发明的一方面，提出了一种增强反Q值补偿稳定性的地震资料处理方法，包括：步骤1：获取地震数据、反射系数和Q值之间的关系表达式；步骤2：获得Lp范数约束下的反演目标函数；步骤3：获得所述Lp范数约束下的反演目标函数的Q值；步骤4：基于所述关系表达式和所述Lp范数约束下的反演目标函数的Q值，获得反射系数。
优选的，所述步骤1包括：子步骤101：基于地震数据建立波动方程，求取所述波动方程的平面波解析解；子步骤102：将所述平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；子步骤103：对所述衰减介质下的平面波解析解进行离散化，获得所述地震数据、反射系数和Q值之间的关系表达式。
优选的，所述衰减介质下的平面波解析解为：
Figure BDA0001769418430000021
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000022
为地震子波，t为时间，
Figure BDA0001769418430000023
为振幅衰减项，
Figure BDA0001769418430000024
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000031
优选的，所述地震数据、Q值和反射系数之间的关系表达式为：
s＝FHWA(Q)r (2)
其中，
Figure BDA0001769418430000032
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000033
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000034
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
优选的，所述步骤2包括：子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于所述Q值建立时变子波库；子步骤202：基于所述时变子波库，建立LP范数约束下的反演目标函数。
优选的，所述LP范数中的P的数值范围为0-1。
优选的，所述Lp范数约束下的反演目标函数为：
Figure BDA0001769418430000035
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000041
为傅里叶逆变换矩阵，
Figure BDA0001769418430000042
为衰减矩阵，
Figure BDA0001769418430000043
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000044
Q为介质的品质因子，W为初始子波矩阵。
根据本发明的另一方面，提出了一种增强反Q值补偿稳定性的地震资料处理系统，该系统包括：存储器，存储有计算机可执行指令；处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：步骤1：获取地震数据、反射系数和Q值之间的关系表达式；步骤2：获得Lp范数约束下的反演目标函数；步骤3：获得所述Lp范数约束下的反演目标函数的Q值；步骤4：基于所述关系表达式和所述Lp范数约束下的反演目标函数的Q值，获得反射系数。
优选的，所述步骤1包括：子步骤101：基于地震数据建立波动方程，求取所述波动方程的平面波解析解；子步骤102：将所述平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；子步骤103：对所述衰减介质下的平面波解析解进行离散化，获得所述地震数据、反射系数和Q值之间的关系表达式。
优选的，所述步骤2包括：子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于所述Q值建立时变子波库；子步骤202：基于所述时变子波库，建立LP范数约束下的反演目标函数。
本发明的有益效果在于：本发明通过获取地震数据、反射系数和Q值之间的关系表达式以及Lp范数约束下的反演目标函数，计算获取Lp范数约束下的反演目标函数的Q值，基于地震数据、反射系数和Q值之间的关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数，增加Lp范数约束的反演方法，较大地提高了地震信号的分辨率和稳定性好。
本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施例中将是显而易见的，或者将在并入本文中的附图和随后的具体实施例中进行详细陈述，这些附图和具体实施例共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的一个实施例的一种增强反Q值补偿稳定性的地震资料处理方法的流程图。
图2a示出了地震数据的吸收衰减记录图。
图2b示出了常规反Q技术的地震数据Q补偿效果图。
图2c示出了阻尼反Q技术的地震数据Q补偿效果图。
图2d示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Q补偿高分辨率处理效果图。
图2e示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的子波整形后的效果图。
图3a示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Lp范数p为1约束下的多次统计图。
图3b示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Lp范数p为0.5约束下的多次统计图。
图3c示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Lp范数p为0.1约束下的多次统计图。
图4a示出了实际地震数据处理效果图。
图4b示出了常规Q补偿高分辨处理的效果图。
图4c示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的处理效果图。
具体实施方式
下面将更详细地描述本发明的优选实施方式。虽然以下描述了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
根据本发明的基于增强反Q值补偿稳定性的地震资料处理方法，包括：步骤1：获取地震数据、反射系数和Q值之间的关系表达式；步骤2：获得Lp范数约束下的反演目标函数；步骤3：获得Lp范数约束下的反演目标函数的Q值；步骤4：基于关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数。
具体的，通过获取地震数据、反射系数和Q值之间的关系表达式以及Lp范数约束下的反演目标函数，计算获取Lp范数约束下的反演目标函数的Q值，基于地震数据、反射系数和Q值之间的关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数。
根据示例性的实施方式基于增强反Q值补偿稳定性的地震资料处理方法增加Lp范数约束的反演方法，较大地提高了地震信号的分辨率和稳定性好。
作为优选方案，步骤1包括：子步骤101：基于地震数据建立波动方程，求取波动方程的平面波解析解；子步骤102：将平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；子步骤103：对衰减介质下的平面波解析解进行离散化，获得地震数据、反射系数和Q值之间的关系表达式。
作为优选方案，衰减介质下的平面波解析解为：
Figure BDA0001769418430000071
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000072
为地震子波，t为时间，
Figure BDA0001769418430000073
为振幅衰减项，
Figure BDA0001769418430000074
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000075
作为优选方案，地震数据、Q值和反射系数之间的关系表达式为：
s＝FHWA(Q)r (2)
其中，
Figure BDA0001769418430000076
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000077
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000081
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
具体的，根据地震数据建立波动方程，求取取波动方程的平面波解析解，再求取衰减介质下的平面波解析解，最后对衰减介质下的平面波解析解进行离散化，获得地震数据、反射系数和Q值之间的关系表达式。主要计算流程如下：
(1)从经典的声波波动方程出发，获取其平面波解析解如下
Figure BDA0001769418430000082
其中，ω为角频率，t为时间，
Figure BDA0001769418430000083
为简谐波，t′＝h/v为传播时间，h为深度，v为速度，r(t′)为反射系数。该解析解可看成不同频率简谐波的一个线性加权叠加，权系数为反射系数。
(2)复速度替换实速度，获取衰减介质下的平面波解析解如下
Figure BDA0001769418430000084
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000085
为地震子波，t为时间，
Figure BDA0001769418430000086
为振幅衰减项，
Figure BDA0001769418430000087
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000088
(3)衰减介质数据解析表达离散化，获取数据与Q值和反射系数之间的数学关系如下：
Figure BDA0001769418430000091
或s＝FHWA(Q)r，
其中，
Figure BDA0001769418430000092
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000093
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000094
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
作为优选方案，步骤2包括：子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于Q值建立时变子波库；子步骤202：基于时变子波库，建立LP范数约束下的反演目标函数。
具体的，通过将地层正Q滤波算子结合到传统稳态褶积模型中，推导非稳态数据时变褶积模型，建立时变子波与Q的关系模型，并基于Q值建立时变子波库，基于时变子波库，建立LP范数约束下的反演目标函数。
在众多反褶积方法中，基于褶积模型的地震数据匹配项始终贯穿于整个反演求解过程，随后的工作基本围绕在如何减少反演多解性问题展开。其中对于反射系数先验信息的描述通过统计学获得并附加于求解过程中，根据统计学原理，相应的数组范数能够反映其作为随机变量的分布情况。其中L∞表征数组满足均匀分布，L2范数对应Gaussian分布，L1范数对应指数分布以及L0范数对应稀疏分布。基于反射系数和噪声序列的分布关系，Debeye和van Riel将反褶积方法统称为Lp范数反褶积。其中包括三种组合，分别为：L1范数反褶积(反射系数和噪声的L1范数最小)；混合范数反褶积(反射系数的L1范数和噪声的L2范数最小)；L2范数反褶积(反射系数和噪声的l2范数最小)。但是由于地震子波的带限性质、时变特征及求解维度的庞大，即便加入关于解的先验信息，最终仍然无法摆脱局部最小值的干扰而不能达到理想的结果。
L2范数反褶积所得到的解是光滑的，且随着深度的增加分辨率是降低的。而很大信号本身就具备稀疏特性(即只包含很少的非零元素)，因此为了克服二范数分辨率较低的缺陷，以及深层分辨的不一致性，本发明求反演过程中的Q值。随着图像重复原及分类等压缩感知领域内技术的发展，基于非凸Lp范数的最小化方法具有稳定性好，收敛快的特点，因此得到广泛应用。
作为优选方案，LP范数中的P的数值范围为0-1。
作为优选方案，Lp范数约束下的反演目标函数为：
Figure BDA0001769418430000101
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000102
为傅里叶逆变换矩阵，
Figure BDA0001769418430000111
为衰减矩阵，
Figure BDA0001769418430000112
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000113
Q为介质的品质因子，W为初始子波矩阵。
具体的，对于叠后非稳态数据，Lp范数约束下的反演目标函数可以表达为：
Figure BDA0001769418430000114
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000115
为傅里叶逆变换矩阵，
Figure BDA0001769418430000116
为衰减矩阵，
Figure BDA0001769418430000117
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000118
Q为介质的品质因子，W为初始子波矩阵。
根据本发明的一种增强反Q值补偿稳定性的地震资料处理系统，该系统包括：存储器，存储有计算机可执行指令；处理器，处理器运行存储器中的计算机可执行指令，执行以下步骤：步骤1：获取地震数据、反射系数和Q值之间的关系表达式；步骤2：获得Lp范数约束下的反演目标函数；步骤3：获得Lp范数约束下的反演目标函数的Q值；步骤4：基于关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数。
具体的，通过获取地震数据、反射系数和Q值之间的关系表达式以及Lp范数约束下的反演目标函数，计算获取Lp范数约束下的反演目标函数的Q值，基于地震数据、反射系数和Q值之间的关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数。
根据示例性的实施方式基于增强反Q值补偿稳定性的地震资料处理系统增加Lp范数约束的反演方法，较大地提高了地震信号的分辨率和稳定性好。
作为优选方案，步骤1包括：子步骤101：基于地震数据建立波动方程，求取波动方程的平面波解析解；子步骤102：将平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；子步骤103：对衰减介质下的平面波解析解进行离散化，获得地震数据、反射系数和Q值之间的关系表达式。
作为优选方案，衰减介质下的平面波解析解为：
Figure BDA0001769418430000121
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000122
为地震子波，t为时间，
Figure BDA0001769418430000123
为振幅衰减项，
Figure BDA0001769418430000124
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000131
作为优选方案，地震数据、Q值和反射系数之间的关系表达式为：
s＝FHWA(Q)r (2)
其中，
Figure BDA0001769418430000132
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000133
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000134
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
具体的，根据地震数据建立波动方程，求取取波动方程的平面波解析解，再求取衰减介质下的平面波解析解，最后对衰减介质下的平面波解析解进行离散化，获得地震数据、反射系数和Q值之间的关系表达式。主要计算流程如下：
(1)从经典的声波波动方程出发，获取其平面波解析解如下
Figure BDA0001769418430000135
其中，ω为角频率，t为时间，
Figure BDA0001769418430000136
为简谐波，t′＝h/v为传播时间，h为深度，v为速度，r(t′)为反射系数。该解析解可看成不同频率简谐波的一个线性加权叠加，权系数为反射系数。
(2)复速度替换实速度，获取衰减介质下的平面波解析解如下
Figure BDA0001769418430000141
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000142
为地震子波，t为时间，
Figure BDA0001769418430000143
为振幅衰减项，
Figure BDA0001769418430000144
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000145
(3)衰减介质数据解析表达离散化，获取数据与Q值和反射系数之间的数学关系如下：
Figure BDA0001769418430000146
或s＝FHWA(Q)r，
其中，
Figure BDA0001769418430000147
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000148
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000151
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
作为优选方案，步骤2包括：子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于Q值建立时变子波库；子步骤202：基于时变子波库，建立LP范数约束下的反演目标函数。
具体的，通过将地层正Q滤波算子结合到传统稳态褶积模型中，推导非稳态数据时变褶积模型，建立时变子波与Q的关系模型，并基于Q值建立时变子波库，基于时变子波库，建立LP范数约束下的反演目标函数。
在众多反褶积方法中，基于褶积模型的地震数据匹配项始终贯穿于整个反演求解过程，随后的工作基本围绕在如何减少反演多解性问题展开。其中对于反射系数先验信息的描述通过统计学获得并附加于求解过程中，根据统计学原理，相应的数组范数能够反映其作为随机变量的分布情况。其中L∞表征数组满足均匀分布，L2范数对应Gaussian分布，L1范数对应指数分布以及L0范数对应稀疏分布。基于反射系数和噪声序列的分布关系，Debeye和van Riel将反褶积方法统称为Lp范数反褶积。其中包括三种组合，分别为：L1范数反褶积(反射系数和噪声的L1范数最小)；混合范数反褶积(反射系数的L1范数和噪声的L2范数最小)；L2范数反褶积(反射系数和噪声的l2范数最小)。但是由于地震子波的带限性质、时变特征及求解维度的庞大，即便加入关于解的先验信息，最终仍然无法摆脱局部最小值的干扰而不能达到理想的结果。
L2范数反褶积所得到的解是光滑的，且随着深度的增加分辨率是降低的。而很大信号本身就具备稀疏特性(即只包含很少的非零元素)，因此为了克服二范数分辨率较低的缺陷，以及深层分辨的不一致性，本发明求反演过程中的Q值。随着图像重复原及分类等压缩感知领域内技术的发展，基于非凸Lp范数的最小化方法具有稳定性好，收敛快的特点，因此得到广泛应用。
作为优选方案，LP范数中的P的数值范围为0-1。
作为优选方案，Lp范数约束下的反演目标函数为：
Figure BDA0001769418430000161
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000162
为傅里叶逆变换矩阵，
Figure BDA0001769418430000163
为衰减矩阵，
Figure BDA0001769418430000164
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000165
Q为介质的品质因子，W为初始子波矩阵。
具体的，对于叠后非稳态数据，Lp范数约束下的反演目标函数可以表达为：
Figure BDA0001769418430000166
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000171
为傅里叶逆变换矩阵，
Figure BDA0001769418430000172
为衰减矩阵，
Figure BDA0001769418430000173
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000174
Q为介质的品质因子，W为初始子波矩阵。
实施例
图1示出了根据本发明的一个实施例的一种增强反Q值补偿稳定性的地震资料处理方法的流程图。
如图1所示，基于增强反Q值补偿稳定性的地震资料处理方法，包括：
步骤1：获取地震数据、反射系数和Q值之间的关系表达式；
步骤1包括子步骤101-子步骤103：
子步骤101：基于地震数据建立波动方程，求取波动方程的平面波解析解；
从经典的声波波动方程出发，获取其平面波解析解如下
Figure BDA0001769418430000175
其中，ω角频率，t为时间，
Figure BDA0001769418430000176
为简谐波，t′＝h/v为传播时间，h为深度，v为速度，r(t′)为反射系数。该解析解可看成不同频率简谐波的一个线性加权叠加，权系数为反射系数。
子步骤102：将平面波将解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；
复速度替换实速度，获取衰减介质下的平面波解析解如下
Figure BDA0001769418430000181
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure BDA0001769418430000182
为地震子波，t为时间，
Figure BDA0001769418430000183
为振幅衰减项，
Figure BDA0001769418430000184
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure BDA0001769418430000185
子步骤103：对衰减介质下的平面波解析解进行离散化，获得地震数据、反射系数和Q值之间的关系表达式；
衰减介质数据解析表达离散化，获取数据与Q值和反射系数之间的数学联系如下
Figure BDA0001769418430000186
或s＝FHWA(Q)r，
其中，
Figure BDA0001769418430000191
为傅里叶算子，H为复数，t为时间，f为频率，
Figure BDA0001769418430000192
为初始子波矩阵，w(f1)为初始子波，
Figure BDA0001769418430000193
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
步骤2：获得Lp范数约束下的反演目标函数；
步骤2包括子步骤201-子步骤202：
子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于Q值建立时变子波库；
子步骤202：基于时变子波库，建立LP范数约束下的反演目标函数；
在众多反褶积方法中，基于褶积模型的地震数据匹配项始终贯穿于整个反演求解过程，随后的工作基本围绕在如何减少反演多解性问题展开。其中对于反射系数先验信息的描述通过统计学获得并附加于求解过程中，根据统计学原理，相应的数组范数能够反映其作为随机变量的分布情况。其中L∞表征数组满足均匀分布，L2范数对应Gaussian分布，L1范数对应指数分布以及L0范数对应稀疏分布。基于反射系数和噪声序列的分布关系，Debeye和van Riel将反褶积方法统称为Lp范数反褶积。其中包括三种组合，分别为：L1范数反褶积(反射系数和噪声的L1范数最小)；混合范数反褶积(反射系数的L1范数和噪声的L2范数最小)；L2范数反褶积(反射系数和噪声的l2范数最小)。但是由于地震子波的带限性质、时变特征及求解维度的庞大，即便加入关于解的先验信息，最终仍然无法摆脱局部最小值的干扰而不能达到理想的结果。
L2范数反褶积所得到的解是光滑的，且随着深度的增加分辨率是降低的。而很大信号本身就具备稀疏特性(即只包含很少的非零元素)，因此为了克服二范数分辨率较低的缺陷，以及深层分辨的不一致性，本发明求反演过程中的Q值。随着图像重复原及分类等压缩感知领域内技术的发展，基于非凸Lp范数的最小化方法具有稳定性好，收敛快的特点，因此得到广泛应用。
其中，LP范数中的P的数值范围为0-1。
其中，Lp范数约束下的反演目标函数为：
Figure BDA0001769418430000201
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure BDA0001769418430000202
为傅里叶逆变换矩阵，
Figure BDA0001769418430000203
为衰减矩阵，
Figure BDA0001769418430000204
为衰减值，t0为初始时间，ω0为参考角频率，ωr为r点处的角频率，γ为衰减项，
Figure BDA0001769418430000205
Q为介质的品质因子，W为初始子波矩阵。
步骤3：获得Lp范数约束下的反演目标函数的Q值；
通过已知参数和参数的初始值设置，对公式(3)求取获得反演目标函数最小值对应的位置，将反演目标函数最小值对应的位置迭代公式(2)中，对公式(2)求取Q值，这个Q值就是最优的Q值，即Lp范数约束下的反演目标函数的Q值。
步骤4：基于关系表达式和Lp范数约束下的反演目标函数的Q值，获得反射系数。
图2a、图2b和图2c分别示出地震数据的吸收衰减记录图、常规反Q技术的地震数据Q补偿效果图、阻尼反Q技术的地震数据Q补偿效果图。图2d和图2e分别示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Q补偿高分辨率处理效果图和子波整形后的效果图。
如图2a、图2b、图2c、图2d和图2e所示，常规反Q技术具有很强的不稳定性，阻尼反Q技术对深层信号Q补偿不足，基于本发明的增强反Q值补偿稳定性的地震资料处理方法能稳定地对浅、中深层信号进行Q补偿，具有较好的效果和较高的分辨率。在图2a-图2e，横坐标表示道数，纵坐标表示时间(ms)。
图3a、图3b和图3c分别示出了根据本发明的一个实施例的增强反Q值补偿稳定性的地震资料处理方法的Lp范数p为1、p为0.5和p为0.1约束下的多次统计图。
如图3a、图3b和图3c所示，Lp范数p为1时，在Q为50的地方有凸性，Lp范数p为0.5时，在Q为50的地方有凸性加强了，Lp范数p为0.1时，在Q为50的地方有凸性进一步加强，为Q值的稳定性提供保障。在图3a-图3c，横坐标表示道数，纵坐标表示时间(ms)。
图4a示出了常规常规Q补偿高分辨处理的效果图。图4b和图4c分别常规常规Q补偿高分辨处理的效果图和增强反Q值补偿稳定性的地震资料处理方法的处理效果图。在图4a-图4c，横坐标表示道数，纵坐标表示时间(s)。
如图4a、图4b和图4c所示，图4a是某工区原始地震数据，分辨率较低，吸收衰减严重，图4b是常规常规Q补偿高分辨处理后的结果，可以看到，效果不是太理想，部分地方同相轴依然粘在一起，图4c是Lp范数约束下的稳定Q补偿高分辨处理效果，分辨率得到极大地提高，且稳定性也比较好。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (7)

1.一种增强反Q值补偿稳定性的地震资料处理方法，其特征在于，包括：
步骤1：获取地震数据、反射系数和Q值之间的关系表达式；
步骤2：获得Lp范数约束下的反演目标函数；
步骤3：获得所述Lp范数约束下的反演目标函数的Q值；
步骤4：基于所述关系表达式和所述Lp范数约束下的反演目标函数的Q值，获得反射系数；
所述步骤2包括：
子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于所述Q值建立时变子波库；
子步骤202：基于所述时变子波库，建立LP范数约束下的反演目标函数；
所述Lp范数约束下的反演目标函数为：
Figure FDA0003067277770000011
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure FDA0003067277770000012
为傅里叶逆变换矩阵，
Figure FDA0003067277770000021
减矩阵，
Figure FDA0003067277770000022
为衰减值，t0为初始时间，ω0为参考角频率，ωd为d点处的角频率，γ为衰减项，
Figure FDA0003067277770000023
Q为介质的品质因子，W为初始子波矩阵。
2.根据权利要求1所述的地震资料处理方法，其特征在于，所述步骤1包括：
子步骤101：基于地震数据建立波动方程，求取所述波动方程的平面波解析解；
子步骤102：将所述平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；
子步骤103：对所述衰减介质下的平面波解析解进行离散化，获得所述地震数据、反射系数和Q值之间的关系表达式。
3.根据权利要求2所述的地震资料处理方法，其特征在于，所述衰减介质下的平面波解析解为：
Figure FDA0003067277770000024
其中，ω为角频率，r(t′)为反射系数，t′为传播时间，
Figure FDA0003067277770000025
为地震子波，t为时间，
Figure FDA0003067277770000031
为振幅衰减项，
Figure FDA0003067277770000032
为相位校正项或子波整形项，ω0为参考角频率，Q(t′)为介质的品质因子，γ为衰减项，
Figure FDA0003067277770000033
4.根据权利要求1所述的地震资料处理方法，其特征在于，所述地震数据、Q值和反射系数之间的关系表达式为：
s＝FHWA(Q)r (2)
其中，
Figure FDA0003067277770000034
为傅里叶算子，H为复数，t为时间，f为频率，
Figure FDA0003067277770000035
为初始子波矩阵，w(f1)为初始子波，
Figure FDA0003067277770000036
为衰减矩阵，α为衰减值，s为地震数据，r为反射系数。
5.根据权利要求1所述的地震资料处理方法，其特征在于，所述LP范数中的P的数值范围为0-1。
6.一种增强反Q值补偿稳定性的地震资料处理系统，其特征在于，该系统包括：
存储器，存储有计算机可执行指令；
处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：
步骤1：获取地震数据、反射系数和Q值之间的关系表达式；
步骤2：获得Lp范数约束下的反演目标函数；
步骤3：获得所述Lp范数约束下的反演目标函数的Q值；
步骤4：基于所述关系表达式和所述Lp范数约束下的反演目标函数的Q值，获得反射系数；
所述步骤2包括：
子步骤201：在稳态褶积模型中结合地层正Q滤波算子，推导非稳态数据时变褶积模型，建立时变子波与Q值的关系模型，并基于所述Q值建立时变子波库；
子步骤202：基于所述时变子波库，建立LP范数约束下的反演目标函数；
所述Lp范数约束下的反演目标函数为：
Figure FDA0003067277770000041
其中，s为地震数据，r为反射系数，p为范数，G＝F-1WA为时变子波矩阵，
Figure FDA0003067277770000042
为傅里叶逆变换矩阵，
Figure FDA0003067277770000051
为衰减矩阵，
Figure FDA0003067277770000052
为衰减值，t0为初始时间，ω0为参考角频率，ωd为d点处的角频率，γ为衰减项，
Figure FDA0003067277770000053
Q为介质的品质因子，W为初始子波矩阵。
7.根据权利要求6所述的地震资料处理系统，其特征在于，所述步骤1包括：
子步骤101：基于地震数据建立波动方程，求取所述波动方程的平面波解析解；
子步骤102：将所述平面波解析解中的实速度替换为复速度，获取衰减介质下的平面波解析解；
子步骤103：对所述衰减介质下的平面波解析解进行离散化，获得所述地震数据、反射系数和Q值之间的关系表达式。
CN201810942502.4A 2018-08-17 2018-08-17 一种增强反q值补偿稳定性的地震资料处理方法及系统 Active CN110837119B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810942502.4A CN110837119B (zh) 2018-08-17 2018-08-17 一种增强反q值补偿稳定性的地震资料处理方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810942502.4A CN110837119B (zh) 2018-08-17 2018-08-17 一种增强反q值补偿稳定性的地震资料处理方法及系统

Publications (2)

Publication Number Publication Date
CN110837119A CN110837119A (zh) 2020-02-25
CN110837119B true CN110837119B (zh) 2021-09-17

Family

ID=69573818

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810942502.4A Active CN110837119B (zh) 2018-08-17 2018-08-17 一种增强反q值补偿稳定性的地震资料处理方法及系统

Country Status (1)

Country Link
CN (1) CN110837119B (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2691796A2 (en) * 2011-03-29 2014-02-05 Total S.A. A process for characterising the evolution of an oil or gas reservoir over time
CN104007467A (zh) * 2014-04-16 2014-08-27 孙赞东 一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法
WO2016065356A1 (en) * 2014-10-24 2016-04-28 Ion Geophysical Corporation Methods for seismic inversion and related seismic data processing
CN107356965A (zh) * 2017-07-20 2017-11-17 中国石油化工股份有限公司 基于加权叠加噪音压制策略的反射系数反演储层预测方法
CN108037531A (zh) * 2017-11-24 2018-05-15 电子科技大学 一种基于广义全变分正则化的地震反演方法及系统

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7072767B2 (en) * 2003-04-01 2006-07-04 Conocophillips Company Simultaneous inversion for source wavelet and AVO parameters from prestack seismic data
AU2011337162B2 (en) * 2010-12-01 2014-11-06 Exxonmobil Upstream Research Company Primary estimation on OBC data and deep tow streamer data

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2691796A2 (en) * 2011-03-29 2014-02-05 Total S.A. A process for characterising the evolution of an oil or gas reservoir over time
CN104007467A (zh) * 2014-04-16 2014-08-27 孙赞东 一种基于混合范数正则化的叠前三参数反演实现的储层与流体预测方法
WO2016065356A1 (en) * 2014-10-24 2016-04-28 Ion Geophysical Corporation Methods for seismic inversion and related seismic data processing
CN107356965A (zh) * 2017-07-20 2017-11-17 中国石油化工股份有限公司 基于加权叠加噪音压制策略的反射系数反演储层预测方法
CN108037531A (zh) * 2017-11-24 2018-05-15 电子科技大学 一种基于广义全变分正则化的地震反演方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
L1-L2范数联合约束稀疏脉冲反演的应用;王宇;《地球科学-中国地质大学学报》;20090930;第34卷(第5期);第835-840页 *
Nonstationary sparsity-constrained seismic deconvolution;Sun Xue-Kai et al.;《APPLIED GEOPHYSICS》;20141231;第11卷(第4期);第459-467页 *
基于时变子波的品质因子估计;冯玮等;《石油地球物理勘探》;20180228;第53卷(第1期);第136-146页 *

Also Published As

Publication number Publication date
CN110837119A (zh) 2020-02-25

Similar Documents

Publication Publication Date Title
CN112835103B (zh) 自适应去鬼波与宽频准零相位反褶积联合处理方法及系统
CN104808245A (zh) 道集优化处理方法及其装置
CN110646841B (zh) 时变稀疏反褶积方法及系统
Banjade et al. Earthquake accelerogram denoising by wavelet-based variational mode decomposition
CN110687597B (zh) 一种基于联合字典的波阻抗反演方法
CN105676292B (zh) 一种基于二维曲波变换的三维地震数据去噪方法
CN110837119B (zh) 一种增强反q值补偿稳定性的地震资料处理方法及系统
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及系统
CN110749923A (zh) 一种基于范数方程提高分辨率的反褶积方法
CN109212609A (zh) 基于波动方程延拓的近地表噪音压制方法
CN109459788B (zh) 地层品质因子计算方法及系统
CN111352157B (zh) 一种横波静校正方法及系统
CN110794455B (zh) 一种地震波传播能量衰减补偿方法
CN113466936B (zh) 一种断层阴影区crp道集获取方法、装置、设备及存储介质
CN114114422B (zh) 基于方向性多尺度分解的叠前地震数据噪声消除方法
CN112764095B (zh) 基于广义能量比的vsp数据q值计算方法及系统
CN117233839B (zh) 地震数据大地吸收衰减三维空间质控方法、系统以及设备
CN109085649B (zh) 一种基于小波变换优化的地震资料去噪方法
CN112099079B (zh) 自适应分频串联反射率反演方法及系统
CN109977350B (zh) 一种超声涡流信号自适应径向基函数神经网络对消方法
CN112415587A (zh) 储层地震波衰减特性分析方法及储层反射系数的反演方法
CN109242781B (zh) 保地质边界的地震图像去噪方法及计算机可读存储介质
CN114185094A (zh) Rms-svd多次波压制方法、装置、电子设备及介质
CN117930358A (zh) 一种岩石超声测试信号去噪方法、装置及电子设备
CN114740530A (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