CN111208561B - 基于时变子波与曲波变换约束的地震声波阻抗反演方法 - Google Patents

Abstract

本发明涉及一种基于时变子波与曲波变换约束的地震声波阻抗反演方法，属于地球物理勘探技术领域，具体包括以下步骤：对二维零偏移距地震数据进行改进广义S变换，并利用谱模拟技术由时频谱提取时变子波；根据反射系数的稀疏性，基于时变子波进行L1范数约束的反射系数反演，得到高精度的反射系数数据；根据连续地层模型，以反射系数反演结果为输入数据，构建并求解基于曲波变换约束的多道声波阻抗反演目标函数，获得横向增强的相对声波阻抗剖面。本发明方法采取稀疏表示与稀疏约束正则化的策略，利用了曲波变换的多尺度性和多方向性，同时充分利用了地震子波的时变性，从而保证反演出的相对声波阻抗结果能够为后续地下地层岩性解释提供可靠的数据基础。

Description

基于时变子波与曲波变换约束的地震声波阻抗反演方法

技术领域

本发明属于地球物理勘探技术领域，涉及一种应用于油气资源勘探的声波阻抗反演方法，特别是关于一种基于时变子波与曲波变换约束的地震声波阻抗反演方法。

背景技术

地震勘探是探测石油、天然气和矿物等自然资源的常规手段，在调查地下构造情况、预测油气圈闭等方面取得了丰硕的成果，而反射系数/声波阻抗的精确反演对于地震地质解释具有重要意义，其反演结果与井资料结合可为解释人员提供准确的岩性分布信息。自上世纪70年代声波阻抗反演出现以来，该项方法技术便受到了地球物理学界的普遍关注。依据反演理论与介质模型的不同，声波阻抗反演可分为基于波动方程和基于褶积模型两大类反演技术。相对前者而言，基于褶积模型的声波阻抗反演技术具有计算量小、易实现的优点，因此，在实际生产中得以广泛应用。

基于褶积模型的声波阻抗反演其传统实现方法是将反射系数作为输入数据，利用递推公式反演获取声波阻抗结果。该方法可以较完整地保留地震数据中的反射特征，将地震反射系数数据映射到地层岩性，从而与录/测井数据进行联合解释，揭示地震记录与地下地层地质特征之间的关系，对储层进行识别描述。然而，这种实现方法对噪声较为敏感，地震数据含噪情况下会使得反演结果产生较大偏差，且起始时刻的声波阻抗值难以准确求取。另外，由于实际地震子波(以下简称子波)受地层吸收衰减等因素的影响会随传播时间及空间变化，对于非平稳地震记录而言，采用时不变子波反演所获得的反射系数数据往往精度较差，进而会使得声波阻抗反演结果不准确。

值得注意的是，由于传统方法采用逐道反演计算模式，若存在地震道缺失现象，也会对反射系数/声波阻抗反演结果产生影响。除以上问题外，逐道反演会使得反演结果的横向连续性被忽略，反射系数/声波阻抗估算值仅与对应道的地震响应有关，这与沉积学原理不符合，即地下介质被层状结构控制，地震记录中的反射波同相轴一般可表现出良好的空间相关性。但是，由于这种传统方法的逐道反演处理模式无法使用阻抗剖面的空间正则化约束，导致反演结果的横向连续性和平滑性不足，会影响最终解释成果的可靠性与精度。

发明内容

本发明要解决的技术问题在于提供一种基于时变子波与曲波变换约束的地震声波阻抗反演方法，本发明利用稀疏表示和稀疏约束正则策略提高地震数据的反演精度，首先采用改进广义S变换对二维零偏移距地震数据进行时频变换，并由时频谱提取时变子波，然后通过L1范数约束的时变反射系数反演获得高精度反射系数剖面，最后基于曲波变换实现多道声波阻抗数据的获取，从而可为后续地震地质解释提供丰富准确的地下岩性分布特征。

本发明采取以下技术方案：

基于时变子波与曲波变换约束的地震声波阻抗反演方法，具体包括以下步骤：

(1)获得地震数据，通过改进广义S变换获取二维零偏移距地震数据的时频特征，并利用谱模拟技术提取时变子波；

a、对时间域非平稳地震记录沿时间方向按下式进行改进广义S变换，得到地震记录相应的时频谱：

其中，MGST(S(τ,f))表示非平稳地震记录s(t)的改进广义S变换时频谱，f表示频率，p(f)和q为改造窗函数的参数，p(f)＝a+bf随频率变化，应用时保持参数q＝1不变，根据不同类型信号和用途调节a和b，使p随频率线性减小，i表示虚数单位，t表示时间，τ代表窗函数的时间位置；

b、采用谱模拟法，即对某时刻τ＝T的地震记录振幅谱|MGST(S(T,f))|进行多项式拟合，得到对应时刻的子波振幅谱|MGST(W(T,f))|，再改变τ的取值，逐点计算子波振幅谱|MGST(W(τ,f))|，实现“时变子波”时频谱的提取；

c、若相位恒为零，则拟合出的“时变子波”各时刻振幅谱经傅里叶反变换后，即是所求零相位时变子波；而当子波为混合相位时，利用希尔伯特变换获得子波最小相位谱，通过Z变换进一步求得相同振幅谱的全部相位谱系列，然后采用子波相位谱扫描技术进行子波的混合相位提取；

(2)根据地层的构造特点，认为反射系数序列具有稀疏性，基于准确提取的时变子波，构建矩阵向量化的多道时变反射系数反演的目标函数，其无约束的拉格朗日形式为：

其中，

表示最优的反射系数反演结果，

表示使其括号内表达式达到最小值时变量

的取值；用符号m＝1,…,M表示地震数据的道数，

表示时间域多道地震数据向量化后构成的列向量，s_m＝[s_1,m,…,s_n,m,…,s_N,m]^T表示第m道地震记录观测值组成的列向量，s_n,m为第m个地震道第n时刻的观测值，T表示矩阵的转置，n＝1,…,N表示时间采样点数，BlockDiag表示分块对角化矩阵，

表示第m道反射系数序列对应的时变子波褶积矩阵，该道地震记录的每个时刻均对应一个长度为L的子波

换而言之，第m个地震道的子波褶积矩阵W_m包含N个随各时刻变化的时变子波

即每个时刻的子波均可不同，w代表子波其各时刻的数值，

表示时间域多道反射系数序列向量化后构成的列向量，r_m＝[r_1,m,…,r_n,m,…,r_N,m]^T表示第m道地震记录对应的反射系数序列，r_n,m为第m个地震道第n时刻的反射系数值；||·||₁为L1范数正则项，拉格朗日乘子λ用于平衡所求反射系数的准确度和稀疏度之间的关系；利用原始对偶对数障碍算法求解式(2)所示目标函数，获得高精度的反射系数数据；

(3)连续地层模型假设弹性参数随深度变化连续可微，即相邻地层之间阻抗差异较小，则第m道反射系数序列r_m与其声波阻抗对数序列ln(z_m)存在线性关系：

r_m＝(1/2)D ln(z_m)，其矩阵-向量乘积的形式为

其中，m＝1,…,M表示地震数据的道数，z_m表示第m道声波阻抗序列，D代表差分矩阵，N代表最大时间采样点；

根据上述连续地层模型中反射系数与声波阻抗二者关系，以时变反射系数反演结果为输入数据，考虑到地震道之间的空间相关性，构建基于曲波变换约束的多道声波阻抗反演目标函数，其可写为基追踪降噪问题：

其中，

表示计算所有矩阵Z，使其后的表达式||C(Z)||₁取得最小值，s.t.表示目标函数(前者)满足约束条件(后者)，Z＝[ln(z₁),…,ln(z_m),…,ln(z_M)]表示多道声波阻抗对数序列按道排列的矩阵，ln(z_m)＝[ln(z_1,m),…,ln(z_N,m)]^T为第m道声波阻抗对数序列组成的列向量，T表示矩阵的转置，C表示曲波变换算子，||C(Z)||₁表示声波阻抗对数值在曲波变换域的L1范数稀疏约束项，||·||_F表示矩阵的Frobenius范数(即矩阵中各项元素的绝对值平方和的开方)，||·||₁为L1范数正则项，D为差分矩阵，R＝[r₁,…,r_m,…,r_M]表示多道反射系数序列按道排列的矩阵，ε为误差参数，其与噪声水平有关；目标函数将曲波变换算子用于L1范数约束反演，即稀疏约束项为||C(Z)||₁，以提升反演解在曲波变换域的稀疏度；利用谱投影梯度算法求解式(3)获得阻抗对数值，取其指数即可求得横向连续性增强的相对声波阻抗数据。

进一步，所述步骤(1)为基于改进广义S变换提取时变子波，当利用谱模拟法由地震记录时频谱拟合“时变子波”时频谱时，假设在时刻τ的子波振幅谱是单峰光滑曲线，构建其频率域数学模型为：

其中，|W(τ,f)|表示子波的振幅谱，f表示频率，k表示常数，i表示阶数，0≤k≤3，2≤I≤7，通过合理设置参数k和I的值，经过多项式拟合求出a_i(i＝0,1,…,I)后将其代入上述模型，即可得到一条光滑的拟合曲线|W(τ,f)|，改变时刻τ的值，按时间顺序估算子波振幅谱。

进一步，所述步骤(2)为L1范数约束的时变反射系数反演，求解式(2)所示目标函数时，参数λ的选取至关重要，在有井资料地区，若反演的反射系数数据与声波测井数据得到的反射界面的标定对比具有良好的对应关系，则表明反演过程使用的参数λ是合理的，反之，则应调整参数λ值重新进行计算匹配；在无井情况下，将反演的反射系数与子波褶积合成地震记录，将其与原始地震数据作差值计算，根据差值剖面中的误差能量大小及随机分布情况定性讨论反演结果的有效性，据此对参数λ进行调整，重复迭代出最佳参数反演结果。

进一步，所述步骤(3)为基于曲波变换约束的多道声波阻抗反演，由于地震资料存在带限问题，即其缺失了低频和高频成分，仅通过求解式(3)所示目标函数无法估算得到非常低频的真实阻抗信息，即所得为相对声波阻抗数据。实际处理时，往往需要根据声波测井资料、地层划分、地震速度分析或地质模型得到声波阻抗的低频背景模型，再对反演结果进行拓频，或者利用先验信息约束反演结果，以拓宽反演结果的频带宽度。

本发明与现有技术相比的有益效果：

本发明由于采取以上技术方案，其具有以下优点：1、本发明方法利用改进广义S变换获取地震记录的时频信息，由获取的时频谱中提取出更符合实际地震波传播规律的时变子波，并将其用于L1范数约束时变反射系数反演，能够有效地提高反射系数反演结果的准确度和稳定性；2、本发明方法以反射系数结果以及连续地层模型为基础，在反演多道声波阻抗时，既利用了曲波变换的多尺度特征提高反演结果的分辨率，还利用了其多方向性特征有效地改善了声波阻抗反演剖面的横向连续性。同时，本发明方法对初始模型的依赖性较低，并且具有很好的抗噪能力。

附图说明

图1为本发明方法所提供的一种基于时变子波与曲波变换约束的地震声波阻抗反演方法的实施流程图；

图2为真实声波阻抗剖面(400-1400ms时段)(a)、采用时变子波(随着时间由400ms至1400ms，零相位雷克子波主频由50Hz逐渐减小到20Hz)合成的对应地震记录剖面(b)，其中黑线标注位置为地震记录第50道；

图3为合成地震剖面第50道地震记录(a)、标准S变换时频谱(b)、广义S变换时频谱(c)、改进广义S变换时频谱(d)、基于改进广义S变换提取的“时变子波”时频谱(e)，其中黑色虚线标注位置分别为580ms和900ms时刻；

图4为合成地震剖面第50道地震记录580ms时刻的子波提取对比图，其中圆圈标注为原始子波(真实子波)，菱形标注为基于标准S变换提取的子波，方形标注为基于广义S变换提取的子波，三角标注为基于改进广义S变换提取的子波；

图5为合成地震剖面第50道地震记录900ms时刻的子波提取对比图，其中圆圈标注为原始子波(真实子波)，菱形标注为基于标准S变换提取的子波，方形标注为基于广义S变换提取的子波，三角标注为基于改进广义S变换提取的子波；

图6为基于时变子波采用L1范数约束反演得到的反射系数剖面；

图7为声波阻抗低频背景模型；

图8为基于常规递推公式反演得到的声波阻抗剖面(a)、本发明方法得到的声波阻抗剖面(b)；

图9为合成地震剖面(图2(b))中加入40％随机噪声后，基于常规递推公式反演得到的声波阻抗剖面(a)、本发明方法反演得到的声波阻抗剖面(b)。

具体实施方式

下面通过实施例来对本发明的技术方案做进一步解释，但本发明的保护范围不受实施例任何形式上的限制。

本发明提出的基于时变子波与曲波变换约束的地震声波阻抗反演方法，利用褶积模型与时变子波褶积矩阵，反演了L1范数约束的反射系数数据，并根据高精度的反射系数反演结果、连续模型假设中反射系数与声波阻抗二者关系，以及曲波变换其近乎最优的稀疏表示优势，构建并求解了地震声波阻抗反演的目标函数，其反演结果的准确性与横向连续性均能够得到明显提高。本发明的基于时变子波与曲波变换约束的地震声波阻抗反演方法主要分为三步，第一步，基于二维零偏移距地震数据，通过改进广义S变换获取其时频特征，并利用谱模拟技术提取时变子波；第二步，利用准确提取的子波构建时变子波褶积矩阵以及相应的L1范数约束目标函数，实现高精度的反射系数反演；第三步，以获得的反射系数结果为输入数据，依据连续地层模型假设，构建基于曲波变换约束的多道声波阻抗反演目标函数，采用谱投影梯度算法进行求解，获得横向增强的相对声波阻抗剖面。

图1为本发明方法的基于时变子波与曲波变换约束的地震声波阻抗反演方法的实施流程图。图2为真实声波阻抗剖面，以及相应的合成地震剖面，采用的子波随着时间变化，具体为雷克子波主频由50Hz逐渐减小到20Hz，剖面记录道数为100，每道样点数为500，采样间隔为2ms。图3～图9为针对图2所示地震数据的声波阻抗反演的相关图件，以下结合附图详细阐述本发明基于时变子波与曲波变换约束的地震声波阻抗反演方法的具体实施过程：

(1)利用改进广义S变换对原始地震数据(见图2(b))逐道进行时频分析，获得对所有频率成分保有较好聚焦性的时频谱。本实施例中对合成地震记录沿时间方向按下式进行改进广义S变换，得到地震记录相应的时频谱(见图3(d))：

其中，MGST(S(τ,f))表示地震记录s(t)的改进广义S变换时频谱，f表示频率，p(f)和q为改造窗函数的参数，p(f)＝a+bf随频率变化，应用时保持参数q＝1不变，根据不同类型信号和用途调节a和b，使p随频率线性减小，i表示虚数单位，s(t)为本实施例合成地震记录，t表示时间，τ代表窗函数的时间位置；

再以本方法获得的地震记录时频谱为基础，采用谱模拟法由逐道时频分析获得的时频谱提取“时变子波”时频谱，即对某时刻τ＝T的地震记录振幅谱|MGST(S(T,f))|进行多项式拟合，得到对应时刻的子波振幅谱|MGST(W(T,f))|，再改变τ的取值，逐点计算子波振幅谱|MGST(W(τ,f))|，实现“时变子波”时频谱的提取，图3(e)为采用谱模拟技术提取的第50道地震记录对应的“时变子波”时频谱；

若原地震记录子波的相位恒为零，则拟合出的“时变子波”各时刻振幅谱经傅里叶反变换后，即是所求零相位时变子波；而当子波为混合相位时，利用希尔伯特变换获得子波最小相位谱，通过Z变换进一步求得相同振幅谱的全部相位谱系列，然后采用子波相位谱扫描技术进行子波的混合相位提取。本实施例中合成的地震记录采用零相位雷克子波，属于第一种情形，故只需对各时刻的子波振幅谱进行傅里叶反变换，即求得零相位时变子波(第50道地震记录580ms提取子波结果如图4中三角标注所示、900ms提取子波结果如图5中三角标注所示)；

(2)利用估计的时变子波建立各单道地震记录的子波褶积矩阵：

本实施例中L表示子波长度，取值为61个样点，子波褶积矩阵包含的时变子波个数N为500，道数m的取值范围为[1,100]，

为第m道提取的时变子波。基于褶积模型构建求解相应的多道时变反射系数反演的目标函数：

其中，

表示最优的反射系数反演结果，

表示使其括号内表达式达到最小值时变量

的取值；用符号m＝1,…,M表示地震数据的道数，

表示时间域多道地震数据向量化后构成的列向量，s_m＝[s_1,m,…,s_n,m,…,s_N,m]^T表示第m道地震记录观测值组成的列向量，T表示矩阵的转置，n＝1,…,N表示时间采样点数，BlockDiag表示分块对角化矩阵，

表示时间域多道反射系数序列向量化后构成的列向量，r_m＝[r_1,m,…,r_n,m,…,r_N,m]^T表示第m道地震记录对应的反射系数序列；||·||₁为L1范数正则项，选取合适参数λ，利用原始对偶对数障碍算法求解上式目标函数，获得高精度的反射系数剖面如图6所示。

求解式(2)所示目标函数时，参数λ的选取至关重要，在有井资料地区，若反演的反射系数数据与声波测井数据得到的反射界面的标定对比具有良好的对应关系，则表明反演过程使用的参数λ是合理的，反之，则应调整参数λ值重新进行计算匹配；在无井情况下，将反演的反射系数与子波褶积合成地震记录，将其与原始地震数据作差值计算，根据差值剖面中的误差能量大小及随机分布情况定性讨论反演结果的有效性，据此对参数λ进行调整，重复迭代出最佳参数反演结果。本实施例是基于合成地震记录进行处理的，没有实际测井资料做比对，故采取第二种策略完成参数λ的优选。

(3)以获得的反射系数数据(见图6)为基础，按照下式进行基于曲波变换约束的多道声波阻抗反演：

其中，

表示计算所有矩阵Z，使其后的表达式||C(Z)||₁取得最小值，s.t.表示目标函数(前者)满足约束条件(后者)，Z＝[ln(z₁),…,ln(z_m),…,ln(z_M)]表示多道声波阻抗对数序列按道排列的矩阵，ln(z_m)＝[ln(z_1,m),…,ln(z_N,m)]^T为第m道声波阻抗对数序列组成的列向量，T表示矩阵的转置，C表示曲波变换算子，||C(Z)||₁表示声波阻抗对数值在曲波变换域的L1范数稀疏约束项，||·||_F表示矩阵的Frobenius范数，||·||₁为L1范数正则项，D为差分矩阵，R＝[r₁,…,r_m,…,r_M]表示多道反射系数序列按道排列的矩阵，ε为误差参数，其与噪声水平有关；利用谱投影梯度算法求解上式获得阻抗对数值，取其指数即可求得横向连续性增强的相对声波阻抗数据。

求解式(3)目标函数时，参数ε取值越大，反演结果的横向连续性越强，而单道声波阻抗曲线的阶梯状变化会变弱，反之则情况相反。选择合适参数ε值，能够得到可以清晰反映层状地层边界的高准确度的相对声波阻抗剖面。

由于地震资料带宽有限，即缺失低频和高频成分，仅通过求解上述目标函数无法直接估算得到真实的绝对声波阻抗数据，即所得为相对声波阻抗数据。实际处理时，往往需要根据声波测井资料、地层划分、地震速度分析或地质模型得到声波阻抗的低频背景模型，再对反演结果进行拓频，或者利用先验信息约束反演结果，以拓宽反演结果的频带宽度。本实施例取真实声波阻抗剖面(见图2(a))的0-8Hz低频成分作为低频背景模型(见图7)，对求解得到的相对声波阻抗数据做8Hz以上高通滤波，将滤波结果与低频背景模型融合，即可获得最终的绝对声波阻抗剖面，如图8(b)所示。

图3所示为第50道地震记录及时频谱图，其中，本方法获得的地震道时频谱(见图3(d))频率聚焦效果明显好于基于标准S变换的时频谱(见图3(b))或基于广义S变换的时频谱(见图3(c))。图4和图5分别为地震剖面(见图2(b))第50道地震记录在580ms和900ms两个不同时刻对应的子波，对比可知，同一位置不同时刻的子波波形具有差异，反映了地震记录的时变特征。相比于基于标准S变换以及广义S变换提取的子波，本发明方法估计的时变子波与原始子波波形更为贴合，这说明本方法提取的时变子波可更好地适应地震记录的非平稳信号特征。图6为基于时变子波的L1范数约束反演得到的反射系数剖面，与输入地震剖面相比，反演的反射系数展布特征清楚，并且具有很高的分辨率，其拓宽了原始地震记录的频带范围。图8为拓频(0-8Hz低频模型见图7)后多道声波阻抗反演结果，图8(b)所示地层声波阻抗呈块状显示，对于地层边缘刻画较为清晰，对地层具有良好的解释能力，且其横向连续性好于传统方法反演结果(见图8(a))。另外，在合成地震剖面(见图2(b))中加入40％随机噪声，基于常规递推公式反演得到的声波阻抗结果(见图9(a))受噪声影响较大，剖面中的层状地层不连续，而本发明方法则具有一定的抗噪能力，所得结果(见图9(b))与真实阻抗模型更为吻合。

模拟实验的地震数据反演处理结果表明，本发明方法能够提高地层参数反演处理的精度，可有效解决岩性组合复杂、横向变化快且地层边界刻画不清的问题，进而提升油气勘探的效率和成功率。

Claims (4)

1.基于时变子波与曲波变换约束的地震声波阻抗反演方法，其特征在于它具体包括以下步骤：

(1)获得地震数据，通过改进广义S变换获取二维零偏移距地震数据的时频特征，并利用谱模拟技术提取时变子波；

a、对时间域非平稳地震记录沿时间方向按下式进行改进广义S变换，得到地震记录相应的时频谱：

其中，MGST(S(τ,f))表示非平稳地震记录s(t)的改进广义S变换时频谱，f表示频率，p(f)和q为改造窗函数的参数，p(f)＝a+bf随频率变化，应用时保持参数q＝1不变，根据不同类型信号和用途调节a和b，使p随频率线性减小，i表示虚数单位，t表示时间，τ代表窗函数的时间位置；

b、采用谱模拟法，即对某时刻τ＝T的地震记录振幅谱|MGST(S(T,f))|进行多项式拟合，得到对应时刻的子波振幅谱|MGST(W(T,f))|，再改变τ的取值，逐点计算子波振幅谱|MGST(W(τ,f))|，实现“时变子波”时频谱的提取；

c、若相位恒为零，则拟合出的“时变子波”各时刻振幅谱经傅里叶反变换后，即是所求零相位时变子波；而当子波为混合相位时，利用希尔伯特变换获得子波最小相位谱，通过Z变换进一步求得相同振幅谱的全部相位谱系列，然后采用子波相位谱扫描技术进行子波的混合相位提取；

(2)根据地层的构造特点，认为反射系数序列具有稀疏性，基于准确提取的时变子波，构建矩阵向量化的多道时变反射系数反演的目标函数，其无约束的拉格朗日形式为：

其中，

表示最优的反射系数反演结果，

表示使其括号内表达式达到最小值时变量

的取值；用符号m＝1,…,M表示地震数据的道数，

表示时间域多道地震数据向量化后构成的列向量，s_m＝[s_1,m,…,s_n,m,…,s_N,m]^T表示第m道地震记录观测值组成的列向量，s_n,m为第m个地震道第n时刻的观测值，T表示矩阵的转置，n＝1,…,N表示时间采样点数，BlockDiag表示分块对角化矩阵，

表示第m道反射系数序列对应的时变子波褶积矩阵，该道地震记录的每个时刻均对应一个长度为L的子波

即第m个地震道的子波褶积矩阵W_m包含N个随各时刻变化的不同子波

w代表子波其各时刻的数值，

表示时间域多道反射系数序列向量化后构成的列向量，r_m＝[r_1,m,…,r_n,m,…,r_N,m]^T表示第m道地震记录对应的反射系数序列，r_n,m为第m个地震道第n时刻的反射系数值；||·||₁为L1范数正则项，拉格朗日乘子λ用于平衡所求反射系数的准确度和稀疏度之间的关系；利用原始对偶对数障碍算法求解式(2)所示目标函数，获得高精度的反射系数数据；

(3)连续地层模型假设弹性参数随深度变化连续可微，即相邻地层之间阻抗差异较小，则第m道反射系数序列r_m与其声波阻抗对数序列ln(z_m)存在线性关系：

r_m＝(1/2)Dln(z_m)，其矩阵-向量乘积的形式为

其中，m＝1,…,M表示地震数据的道数，z_m表示第m道声波阻抗序列，D代表差分矩阵，N代表最大时间采样点；

以时变反射系数反演结果为输入数据，考虑到地震道之间的空间相关性，构建基于曲波变换约束的多道声波阻抗反演目标函数，其可写为基追踪降噪问题：

其中，

表示计算所有矩阵Z，使其后的表达式||C(Z)||₁取得最小值，s.t.表示目标函数满足约束条件，Z＝[ln(z₁),…,ln(z_m),…,ln(z_M)]表示多道声波阻抗对数序列按道排列的矩阵，ln(z_m)＝[ln(z_1,m),…,ln(z_N,m)]^T为第m道声波阻抗对数序列组成的列向量，T表示矩阵的转置，C表示曲波变换算子，||C(Z)||₁表示声波阻抗对数值在曲波变换域的L1范数稀疏约束项，||·||_F表示矩阵的Frobenius范数，||·||₁为L1范数正则项，D为差分矩阵，R＝[r₁,…,r_m,…,r_M]表示多道反射系数序列按道排列的矩阵，ε为误差参数，其与噪声水平有关；利用谱投影梯度算法求解式(3)所示目标函数可以获得阻抗对数值，取其指数即可求得横向连续性增强的相对声波阻抗数据。

2.根据权利要求1所述的方法，其特征在于所述步骤(1)为基于改进广义S变换提取时变子波，当利用谱模拟法由地震记录时频谱拟合“时变子波”时频谱时，假设在时刻τ的子波振幅谱是单峰光滑曲线，构建其频率域数学模型为：

其中，|W(τ,f)|表示子波振幅谱，f表示频率，k表示常数，i＝0,1,…,I表示阶数，0≤k≤3，2≤I≤7，通过合理设置参数k和I的值，经过多项式拟合求出a_i后将其代入上述模型，即可得到一条光滑的拟合曲线|W(τ,f)|，a_i表示多项式系数,改变时刻τ的值，按时间顺序估算子波振幅谱。

3.根据权利要求1所述的方法，其特征在于所述步骤(2)为L1范数约束的时变反射系数反演，求解式(2)所示目标函数时，参数λ的选取至关重要，在有井资料地区，若反演的反射系数数据与声波测井数据得到的反射界面的标定对比具有良好的对应关系，则表明反演过程使用的参数λ是合理的，反之，则应调整参数λ值重新进行计算匹配；在无井情况下，将反演的反射系数与子波褶积合成地震记录，将其与原始地震数据作差值计算，根据差值剖面中的误差能量大小及随机分布情况定性地讨论反演结果的有效性，据此对参数λ进行调整，重复迭代出最佳反演结果。

4.根据权利要求1所述的方法，其特征在于所述步骤(3)为基于曲波变换约束的多道声波阻抗反演，由于实际地震资料存在带限问题，即其缺失低频和高频成分，仅通过求解式(3)所示目标函数无法估算得到非常低频的真实阻抗信息，即所得为相对声波阻抗数据；实际处理时，需要根据声波测井资料、地层划分、地震速度分析或地质模型得到声波阻抗的低频背景模型，再对反演结果进行拓频，或者利用目标函数反演解的先验信息约束反演结果，以拓宽反演结果的频带宽度。

Images