基于全变差最小化约束的地震反射系数反演方法
【技术领域】
[0001] 本发明属于地震油气勘探技术领域，具体为一种基于全变差最小化约束共辄梯度 算法的地震反射系数反演方法。
【背景技术】
[0002] 作为地震油气勘探中的重要步骤和环节，地震储层预测一直在油气勘探中扮演着 不可或缺的角色。随着岩性油气藏勘探开发的不断深入，对小于地震波长四分之一的薄储 层(薄层)预测问题日益突出。模型反演，稀疏脉冲反演等常规地震反演方法对薄储层识别 能力有限，直接影响到后续的油气藏评价和井位设计。
[0003] 地震反射系数反演是在频率域实施的，可获得高分辨率的时间域反射系数信息， 实现对小于调谐厚度的薄储层进行有效识别，以达到提高地震数据分辨率和提高储层预测 精度的目的。
[0004] 地震反射系数反演是用有限的频域信息得到全域频域信息，是一个欠定反演问 题，传统共辄梯度方法求解欠定问题效率低且易受异常值影响，会导致反演结果背离实际 地层特性，与地震剖面匹配较差等不利后果，不能正确的反应出真实的地层信息。
[0005] 目前应用于地震反射系数估计的反演算法包括：Charles I .Puryear和John P.Castagna采用的最小二乘共辄梯度法。其基本思想是把共辄性与最小二乘方法相结合， 利用已知点处的梯度构造一组共辄方向，并沿这组方向进行搜素，求出目标函数的极小点。
[0006] Rui Zhang和John P.Castagna采用的基追踪算法。它采用反射系数的范数作为稀 疏性的度量，通过最小化L1范数将反射系数稀疏表示问题定义为一类有约束的极值问题， 进而转化为线性规划问题进行求解。
[0007] 袁三一和王尚旭等采用的粒子群算法和列文伯格-马夸尔特法的联合算法。其首 先采用粒子群算法进行全局寻优，来反演出地层反射系数位置，然后采用列文伯格-马夸尔 特法算法来精准地求取反射系数的数值大小。粒子群算法和列文伯格-马夸尔特法算法相 联合，反演结果更精确，收敛性更快。
[0008] 秦德文等采用的基于随机搜索的蒙特卡罗算法。它是以概率和统计理论方法为基 础的一种计算方法。将所求解的问题同一定的概率模型相联系，用计算机实现统计模拟或 抽样，以获得问题的近似解。蒙特卡罗方法有很强的适应性，该方法的收敛性是指概率意义 下的收敛，因此问题维数的增加不会影响它的收敛速度。
[0009] Kelyn PaolaCastafio.和German Ojeda采用的模拟退火算法和遗传算法。模拟退 火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随 机寻找目标函数的全局最优解。而遗传算法是模拟达尔文生物进化论的自然选择和遗传学 机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。应用 于反射系数估计问题时，遗传算法效果比模拟退火算法稍好。
[0010] 柴新涛，李振春等采用的最小二乘QR分解算法。它是利用Lanezos方法求解最小二 乘问题的一种投影法。由于在求解过程中用到QR因子分解，在对数据误差传递的压制和求 解收敛效率上具有明显的优越性。总体来讲，这些算法都可以不同程度的解决反射系数的 反演问题。
【发明内容】
[0011] 针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的 一个或多个问题。例如，本发明的目的之一在于解决现有技术中存在的反射系数反演效率 低和易受异常值影响的问题，提供一种基于全变差最小化约束的共辄梯度地震反射系数反 演方法。该方法利用地层连续性约束，在频率域提升共辄梯度反演的准确性和效率。通过全 变差最小化约束增加对反射系数的识别能力，抑制共辄梯度方法易受到数据变化扰动以及 反演结果地层跳变的问题，获得高分辨率的时间域反射系数，提升地震数据的分辨率，提高 对小于调谐厚度的薄储层反射系数识别能力。
[0012] 为了实现上述目的，本发明提供了一种基于全变差最小化约束的地震反射系数反 演方法。所述反演方法包括以下步骤：
[0013] A、从叠后地震数据S中提取地震子波w并进行傅立叶变换得到地震子波w的频域表 示W(f)，其中：
[0014] ff(f)=FFT(w) (1)
[0015] B、从所述叠后地震数据S中取一个地震道的一个时窗内的地震数据s进行傅立叶 变换，得到地震数据s的频域表示S(f)，其中：
[0016] S(f)=FFT(s) (2)
[0017] C、根据地震数据形成原理由所述地震数据s的频域表示S(f)和所述地震子波w的 频域表示w(f)得到反射系数的频域表示R(f)，其中：
[0018] R(f)=S(f)/ff(f) (3)
[0019] D、以所述时窗中心为分析点，对所述时窗内的地震数据s的时域反射系数进行傅 立叶变换，得到反射系数的频域表达R(fT，其中：
[0020]
⑷
[0021 ]在等式(4)中，N为由点数表示的分析数据长度，Ti代表反射系数对之间的时间间 隔，j为虚数单位。
[0022] E、在所述步骤C得到的反射系数的频域表示R(f)和所述步骤D得到的反射系数的 频域表达R(fV相等的基础上，根据反射系数奇偶分解原理构建所述时窗内反射系数反演 的目标函数方程。
[0023] F、利用共辄梯度算法求解所述目标函数方程以得到所述时窗内反射系数的奇分 量和偶分量，并重构出所述时窗内的反射系数，其中，在每次迭代过程中对共辄梯度算法所 求得的解进行全变差最小化约束，同时将经全变差最小约束之后的解作为下一次迭代计算 的初始解。
[0024] G、在所述地震道上按步长step滑动所述时窗，得到下一个时窗的地震数据，重复 步骤B至F，直到所述时窗遍历所述地震道的地震数据，完成所述地震道的反演，得到所述地 震道的反射系数。
[0025] H、取下一地震道的地震数据并按照步骤B到G进行处理，直到得到每一个地震道的 反射系数。
[0026] 根据本发明的基于全变差最小化约束的地震反射系数反演方法的一个实施例，所 述步骤E包括:用欧拉公式将所述步骤D得到反射系数的频域表达R(fV中的指数项展开，得 到等式(5):
[0027]
(5)
[0028]根据奇偶分解原理，等式(5)中实部是反射系数对^与^1+1的偶分量，用re(i，N_i +1)来表示，虚部是反射系数对ri与rN-i+i的奇分量，用r〇(i，N-i+l)来表示，得到等式(6):
[0029]
(6)
[0030]将所述等式⑶分析点移至时窗中点，并结合等式(6)建立目标函数方程，得到等式(7):
[0031]
(7 )
[0032] 在等式(7)中，flOT为低频截止频率，fhigh为高频截止频率，~为偶分量权重，α。为奇 分量权重，Re代表着实部，Im代表着虚部，h为反射系数的偶分量，:r。为反射系数的奇分量， △ t为相对于分析点的位移量，tw为时窗win的半窗长。
[0033]将等式(7)在频域和时域离散化并写为奇分量和偶分量的形式，得到下面的等式 (8)和等式(9):
[0034]
[0035]
[0036] 在等式(8)和等式(9)中，是起始频率，取值为所述flOT; fm是截止频率，取值为所 high 〇
[0037] 将所述等式(8)和等式(9)分别写为矩阵形式，得到下面的等式(10)和等式(11):
[0038] be = AeXre (10)
[0039] b〇 = A〇Xr〇 (11)
[0040] 建立目标函数方程：
[0041 ] (12)
[0042] 在等式（10)、（ 11)和（12)中，为偶分量权重，α。为奇分量权重;是频域反射系数 的实部;b。是频域反射系数的虚部;I为偶分量变换矩阵，BP:
[0043]
(13)
[0044] A。为奇分量变换矩阵，即：
[0045]
(.1.4).
[0046] 将所述目标函数方程(12)写为更一般的优化问题，得到目标函数方程(15):
[0047] b=Ax (15)
[0048] 在等式(15)中，X为待求的反射系数。
[0049] 根据本发明的基于全变差最小化约束的地震反射系数反演方法的一个实施例，所 述步骤F包括：
[0050] 将所述目标函数方程（15)的求解问题转化为| |Ax=b| |2优化问题，进一步写成二 次函数形式：
[0051 ]
(16)
[0052] 在等式（16)中，Q=ATA，y = ATb，T表示转置运算。
[0053] 初始化:叉。=0，1'。= 13-厶叉。，〇=||(>)||2，1^=1，1^为迭代次数。
[0054] 循环步骤：
[0055] (1)更新梯度方向：
[0056] gk = Qxk-i-y (17)
[0057] (2)更新共辄方向：
[0058]若k = l，则pk = r。；若k关 1，则
[0060] pk = _gk+akPk-1 (19)[0061] (3)更新第k次迭代的解：
[0059] (18) _2]
(2〇) [0063] (4)对xk进行全变差最小化，
[0064] ,~、 (21)
[0065]其中，D表示整体剖面集，采用对称边界条件，xu表示第i行第j个元素。
[0066]