本发明公开了提高地震分辨率的保构造井控空间反褶积方法和装置。该方法在常规反褶积模型基础上采用沿着空间方向正则化技术，实现多道同时反褶积；将井信息通过初始模型引入到了反演结果中，加以约束和质控，使结果更加准确，避免空间不稳定现象，提高纵向分辨率的同时，保持高分辨率结果的横向连续性。实际应用中，本方法可以有效地提高对薄砂体的识别准确度，对后续的叠后地震属性分析也具有十分重要的实际意义。
1.一种提高地震分辨率的保构造井控空间反褶积方法，该方法包括：
步骤102，给出初始的空间正则化参数序列λi,i＝1,2,...,N，N为原始的多道地震记录中包含的总道数，进入步骤103；
步骤103，采用沿空间方向的正则化，基于当前的空间正则化参数序列对所述原始的多道地震记录进行空间反褶积，得到第二反射系数矩阵，进入步骤104；
步骤104，如果得到的第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，则进入步骤106，否则，进入步骤105；
步骤105，调整当前的空间正则化参数序列，然后回到步骤103；
步骤106，将当前得到的第二反射系数矩阵作为最终的反射系数矩阵，将其中每一道对应的反射系数序列与该道的子波进行褶积，并叠加多道的褶积结果，以得到所需的叠后地震数据体。
2.根据权利要求1所述的方法，其中，在步骤103中，得到第二反射系数矩阵包括：设定目标函数E(m)：
其中，d是对原始的多道地震记录进行叠加后得到的原始叠后地震数据体，G为子波褶积矩阵；m表示待求的第二反射系数矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||2表示2范数；
迭代求解上式得到使得目标函数E(m)小于第一预设阈值的矩阵m，该矩阵m即为在当前的空间正则化参数序列下得到的第二反射系数矩阵。
3.根据权利要求1所述的方法，其中，在步骤104中，通过下列方法判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望：设定目标函数min(m)：
其中，m表示第二反射系数矩阵，S表示第一反射系数矩阵，C表示协方差矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||p表示p范数；
如果min(m)小于第二预设阈值，则判断第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，反之，则判断第一反射系数矩阵和第二反射系数矩阵间的误差不符合预定期望。
4.一种提高地震分辨率的保构造井控空间反褶积装置，该装置包括第一反射系数矩阵确定单元、初始化单元、第二反射系数矩阵确定单元、误差判断单元、空间正则化参数调整单元、褶积叠加单元，其中：第一反射系数矩阵确定单元，用于基于测井曲线得到第一反射系数矩阵；
初始化单元，用于给出初始的空间正则化参数序列λi,i＝1,2,...,N，其中N为原始的多道地震记录中包含的总道数；
第二反射系数矩阵确定单元，采用沿空间方向的正则化，用于基于当前的空间正则化参数序列对所述原始的多道地震记录进行空间反褶积以得到第二反射系数矩阵，开始时从初始化单元接收空间正则化参数序列，在误差判断单元判断误差不符合预定期望的情况下，从空间正则化参数调整单元接收空间正则化参数序列；
误差判断单元，分别从第一反射系数矩阵确定单元和第二反射系数矩阵确定单元接收第一反射系数矩阵和第二反射系数矩阵，用于判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望；
空间正则化参数调整单元，在误差判断单元判断所述误差不符合预定期望的情况下，用于调整当前的空间正则化参数序列；
褶积叠加单元，在误差判断单元判断所述误差符合预定期望的情况下，用于将第二反射系数矩阵确定单元当前得到的第二反射系数矩阵作为最终的反射系数矩阵，并将其中每一道对应的反射系数序列与该道的子波进行褶积，以及叠加多道的褶积结果，以得到所需的叠后地震数据体。
5.根据权利要求4所述的装置，其中，在第二反射系数矩阵确定单元中，得到第二反射系数矩阵包括：设定目标函数E(m)：
其中，d是对原始的多道地震记录进行叠加后得到的原始叠后地震数据体，G为子波褶积矩阵；m表示待求的第二反射系数矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||2表示2范数；
迭代求解上式得到使得目标函数E(m)小于第一预设阈值的矩阵m，该矩阵m即为在当前的空间正则化参数序列下得到的第二反射系数矩阵。
6.根据权利要求4所述的装置，其中，在误差判断单元中，通过下列方法判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望：设定目标函数min(m)：
其中，m表示第二反射系数矩阵，S表示第一反射系数矩阵，C表示协方差矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||p表示p范数；
如果min(m)小于第二预设阈值，则判断第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，反之，则判断第一反射系数矩阵和第二反射系数矩阵间的误差不符合预定期望。
提高地震分辨率的保构造井控空间反褶积方法和装置
技术领域
[0001] 本发明涉及油气地球物理技术领域，更具体地，涉及一种提高地震分辨率的保构造井控空间反褶积方法和一种提高地震分辨率的保构造井控空间反褶积装置。
背景技术
[0002] 叠后地震资料可以看成子波和反射系数的褶积，这导致了子波和资料的同时带限。这意味着子波滤掉了一些地震资料频带外的重要地质信息。高分辨率地震资料处理结果，对于薄储层、微幅构造、小断层等小尺度地质体的识别，开发层系细分、薄油层识别、井间连通性分析，以及建立精细的三维油藏地质模型至关重要。而反褶积是提高地震资料分辨率和拓展频带的一种重要手段。目前的反褶积方法大体上可以分为两类，一类是以预测反褶积为代表的反褶积方法；另一类是不对地震子波的相位做任何假设的反褶积方法，如同态反褶积、最小熵反褶积等。这些常用方法均有其各自的优、缺点。
[0003] 预测反褶积应用于地面地震勘探中压制多次波外，还用于消除海上鸣震作用等方面(胡镠.多道地震记录的预测反褶积算子提取的实现方法及应用.石油仪器.2009，23(5)：55-57)。它是基于线性预测理论的最小平方滤波，即维纳滤波。任何信号本身都存在随机干扰，做信号处理的目的就是要提取包含在随机信号中的确定成分。维纳滤波方法就是这样一种线性滤波器，当信号和随机噪声同时输入该滤波器时能将输出端信号尽可能精确表现出来。它具有两个基本的假设条件，一是地震子波是最小相位，二是反射系数序列是随机过程。但是实际地震资料很难满足这两个假设条件。
[0004] 同态反褶积是一种非线性滤波方法。它应用同态滤波理论，将一种函数映射至另一函数域，使以褶积关系组成的信号先变换为相乘的形式，再把相乘关系组成的信号转换为相加的形式(复塞谱)。它是在复塞谱域内进行滤波运算后，经过反函数转换到原始函数域的广义是、线性滤波系统。但是同态反褶积要求子波以低频为主，反射系数序列以高频为主，因此在复赛谱上能够将两者分离(事实上两者是部分重合的，因而不能将两者完全分离)，从而得到地震子波。王君等(同态反褶积的改进与应用.石油地球物理勘探.2003，38卷：27-30)采用最小熵反褶积MED、L模和PARSIMONY为目标函数，从叠加平均对数谱得到的地震子波和反子波当中自动确定最佳子波和反子波，减少人为因素。然后对地震道做混合相位反褶积，这样就把同态反褶积和最小熵反褶积、L模反褶积结合到一起，相互取长补短。
[0005] 周超等(信噪比谱约束下的井控反褶积方法.岩性油气藏.2010，7：72-76)通过引入测井资料的高频信息，对常规反褶积方法进行改进，实现了井控反褶积；以信噪比谱为质量控制选取约束井控反褶积参数，确保反褶积结果提高分辨率的同时保持高信噪比。该方法可使目的层段地震资料频谱得到一定的拓展，有效带宽内信噪比基本保持不变，并略有提高。与本发明的相似之处在于反褶积过程中均加入了井信息进行质控，但是，当资料的信噪比较低或地震资料包含高波数噪声和误差时，该方法会导致估计的反射系数或阻抗横向不稳定。
[0006] 李志娜等(专利申请号CN104914466A，2015，9月)提出了一种提高地震资料分辨率的方法。该技术为涉及多次波发育地震资料的叠前处理。首先基于SRME和聚焦变换理论推导得出利用多次波的多维子波反褶积模型；然后利用引入整形正则化的非稳态回归自适应匹配滤波方法将聚焦域内由多次波恢复的高分辨率数据分离出来，实现原始数据的高分辨率转换。该技术实现基于波动理论，从空间上加以约束，在提高纵向分辨率的同时可以改善数据的横向分辨率。与本发明的相似之处在于反褶积过程中均从空间上加以约束，但是二者基 于的原理以及实现的方法并不相同，最终要的一点是该专利并未加入井的信息加以约束及质控。
[0007] 总之，目前常规的反褶积方法多应用沿着时间方向添加一个正则化项，并且通过反演方式获得一个稳定和最优的解。但是，这些常规方法大多采用的是一道道处理方式。当资料的信噪比较低或地震资料包含高波数噪声和误差时，单道正则化方法经常会导致估计的反射系数或阻抗横向不稳定。此外，常规的反褶积处理忽略了道与道之间的空间关系，并且未加入井的信息加以约束及质控。
发明内容
[0008] 本发明提出了一种提高地震分辨率的保构造井控空间反褶积方法，本发明还提出了相应的装置。
[0009] 根据本发明的一方面，提出了一种提高地震分辨率的保构造井控空间反褶积方法，该方法包括：
[0010] 步骤101，基于测井曲线得到第一反射系数矩阵；
[0011] 步骤102，给出初始的空间正则化参数序列λi,i＝1,2,...,N，N为原始的多道地震记录中包含的总道数，进入步骤,103；
[0012] 步骤103，采用沿空间方向的正则化，基于当前的空间正则化参数序列对所述多道地震记录进行空间反褶积，得到第二反射系数矩阵，进入步骤104；
[0013] 步骤104，如果得到的第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，则进入步骤106，否则，进入步骤105；
[0014] 步骤105，调整当前的空间正则化参数序列，然后回到步骤103；
[0015] 步骤106，将当前得到的第二反射系数矩阵作为最终的反射系数矩阵，将其中每一道对应的反射系数序列与该道的子波进行褶积，并叠加多道的褶积结果，以得到所需的叠后地震数据体。
[0016] 根据本发明的另一方面，公开了一种提高地震分辨率的保构造井控空间反褶积装置，该装置包括第一反射系数矩阵确定单元、初始化单元、第二反射系数矩阵确定单元、误差判断单元、空间正则化参数调整单元、褶积叠加单元，其中：
[0017] 第一反射系数矩阵确定单元，用于基于测井曲线得到第一反射系数矩阵；
[0018] 初始化单元，用于给出初始的空间正则化参数序列λi,i＝1,2,...,N，其中N为原始的多道地震记录中包含的总道数；
[0019] 第二反射系数矩阵确定单元，采用沿空间方向的正则化，用于基于当前的空间正则化参数序列对所述多道地震记录进行空间反褶积以得到第二反射系数矩阵，开始时从初始化单元接收空间正则化参数序列，在误差判断单元判断所述误差不符合预定期望的情况下，从空间正则化参数调整单元接收空间正则化参数序列；
[0020] 误差判断单元，分别从第一反射系数矩阵确定单元和第二反射系数矩阵确定单元接收第一反射系数矩阵和第二反射系数矩阵，用于判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望；
[0021] 空间正则化参数调整单元，在误差判断单元判断所述误差不符合预定期望的情况下，用于调整当前的空间正则化参数序列；
[0022] 褶积叠加单元，在误差判断单元判断所述误差符合预定期望的情况下，用于将第二反射系数矩阵确定单元当前得到的第二反射系数矩阵作为最终的反射系数矩阵，并将其中每一道对应的反射系数序列与该道的子波进行褶积，以及叠加多道的褶积结果，以得到所需的叠后地震数据体。
[0023] 本发明中，通过保构造井控空间反褶积，解决了单道地震数据正则化方法导致的反射系数或波阻抗横向不稳定的问题；同时将井的信息通过初始模型引入到反演结果中，加以约束和质控，达到提高地震资料分辨率的目的。本发明能够大大提高了地震反褶积信号处理的准确性，具有极高的工业实用价值和推广应用前景。
附图说明
[0024] 通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。
[0025] 图1示出了根据本发明一个实施例的保构造井控空间反褶积方法的流程图。
[0026] 图2示出了根据本发明一个具体应用实例的实现示意图。
[0027] 图3(a)、(b)、(c)示出了某原始含噪地震数据体图、无正则化反褶积叠后地震数据体图、应用本发明得到的叠后地震数据体图。
[0028] 图4示出了应用本发明得到的A井、B井、C井的井控多道反褶积结果与井滤波反射系数对比图。
[0029] 图5(a)、(b)示出了某原始地震剖面与应用本发明得到的反褶积成果剖面对比图。
具体实施方式
[0030] 下面将参照附图更详细地描述本发明的优选实施方式。虽然附图中显示了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
[0031] 实施例1
[0032] 图1示出了根据本发明一个实施例的保构造井控空间反褶积方法的流程图。根据本实施例，该方法包括：
[0033] 步骤101，基于测井曲线得到第一反射系数矩阵，例如，可以在对采集到的测井曲线进行滤波地震信号处理操作后提取井的反射系数序列；
[0034] 步骤102，给出初始的空间正则化参数序列λi,i＝1,2,...,N，N为原始的多道地震记录中包含的总道数，进入步骤,103；
[0035] 步骤103，采用沿空间方向的正则化，基于当前的空间正则化参数序列对所述多道地震记录进行空间反褶积，得到第二反射系数矩阵，进入步骤104；
[0036] 步骤104，如果得到的第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，则进入步骤106，否则，进入步骤105；
[0037] 步骤105，调整当前的空间正则化参数序列，然后回到步骤103；
[0038] 步骤106，将当前得到的第二反射系数矩阵作为最终的反射系数矩阵，将其中每一道对应的反射系数序列与该道的子波进行褶积，并叠加多道的褶积结果，以得到所需的叠后地震数据体。
[0039] 上述实施例至少具有下列有益效果：1.利用空间正则化技术实现多道地震数据同时进行反褶积地震信号处理，改善了传统方法因单道地震数据正则化方法引起的反射系数或阻抗横向不稳定现象；2.使用测井信息进行约束和质控，使得反演结果更准确，更好地保护了地层的边界，像断层、裂缝、尖灭和砂体的边缘等。此外，本领域技术人员可以理解的是，步骤101不必须在步骤102～102之前执行，步骤101可以在步骤104之前的任意时间执行。
[0040] 可选地，在步骤103中，得到第二反射系数矩阵可以包括：
[0041] 设定目标函数E(m)：
[0042]
[0043] 其中，d是对原始的多道地震记录进行叠加后得到的原始叠后地震数据体，G为子波褶积矩阵，m表示待求的第二反射系数矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||2表示2范数，可以被视为实际地震数据与合成地震记录的残差能量， 可以被视为加入的正则化因子，以起到稳定算法的作用；
[0044] 迭代求解上式得到使得目标函数E(m)小于第一预设阈值的矩阵m，该矩阵m即为在当前的空间正则化参数序列下得到的第二反射系数矩阵。
[0045] 可选地，在步骤104中，可以通过下列方法判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望：
[0046] 设定目标函数min(m)：
[0047]
[0048] 其中，m表示第二反射系数矩阵，S表示第一反射系数矩阵，C表示协方差矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||p表示p范数，||C(m-S)||p可以被视为计算出的理论值与井观测数据的差，可以被视为加入的正则化因子，以起到稳定算法的作用；
[0049] 如果min(m)小于第二预设阈值，则判断第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，反之，则判断第一反射系数矩阵和第二反射系数矩阵间的误差不符合预定期望。
[0050] 实施例2
[0051] 本发明还公开了一种提高地震分辨率的保构造井控空间反褶积装置。本实施例中，第一反射系数矩阵确定单元、初始化单元、第二反射系数矩阵确定单元、误差判断单元、空间正则化参数调整单元、褶积叠加单元，其中：
[0052] 第一反射系数矩阵确定单元，用于基于测井曲线得到第一反射系数矩阵；
[0053] 初始化单元，用于给出初始的空间正则化参数序列λi,i＝1,2,...,N，其中N为原始的多道地震记录中包含的总道数；
[0054] 第二反射系数矩阵确定单元，采用沿空间方向的正则化，用于基于当前的空间正则化参数序列对所述多道地震记录进行空间反褶积以得到第二反射系数矩阵，开始时从初始化单元接收空间正则化参数序列，在误差判断单元判断所 述误差不符合预定期望的情况下，从空间正则化参数调整单元接收空间正则化参数序列；
[0055] 误差判断单元，分别从第一反射系数矩阵确定单元和第二反射系数矩阵确定单元接收第一反射系数矩阵和第二反射系数矩阵，用于判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望；
[0056] 空间正则化参数调整单元，在误差判断单元判断所述误差不符合预定期望的情况下，用于调整当前的空间正则化参数序列；
[0057] 褶积叠加单元，在误差判断单元判断所述误差符合预定期望的情况下，用于将第二反射系数矩阵确定单元当前得到的第二反射系数矩阵作为最终的反射系数矩阵，并将其中每一道对应的反射系数序列与该道的子波进行褶积，以及叠加多道的褶积结果，以得到所需的叠后地震数据体。
[0058] 可选地，在第二反射系数矩阵确定单元中，得到第二反射系数矩阵可以包括：
[0059] 设定目标函数E(m)：
[0060]
[0061] 其中，d是对原始的多道地震记录进行叠加后得到的原始叠后地震数据体，G为子波褶积矩阵；m表示待求的第二反射系数矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||2表示2范数；
[0062] 迭代求解上式得到使得目标函数E(m)小于第一预设阈值的矩阵m，该矩阵m即为在当前的空间正则化参数序列下得到的第二反射系数矩阵。
[0063] 可选地，在误差判断单元中，可以通过下列方法判断第一反射系数矩阵和第二反射系数矩阵间的误差是否符合预定期望：
[0064] 设定目标函数min(m)：
[0065]
[0066] 其中，m表示第二反射系数矩阵，S表示第一反射系数矩阵，C表示协方差矩阵，λi为第i道对应的空间正则化参数，Di为第i道对应的表示空间结构的频率-波数域算子，||*||1表示1范数，||*||p表示p范数；
[0067] 如果min(m)小于第二预设阈值，则判断第一反射系数矩阵和第二反射系数矩阵间的误差符合预定期望，反之，则判断第一反射系数矩阵和第二反射系数矩阵间的误差不符合预定期望。
[0068] 应用示例
[0069] 为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
[0070] 图2示出了根据本发明一个具体应用实例的实现示意图。
[0071] 图3(a)、(b)、(c)示出了某原始含噪地震数据体图、无正则化反褶积叠后地震数据体图、应用本发明得到的叠后地震数据体图。
[0072] 图3(a)示出了共计681道地震记录，时间采样间隔是2ms，Amoco模型导出的反射系数与30Hz主频Ricker子波褶积后，加了20％随机噪声作为原始的测试资料。图3(b)示出了无正则化项反褶积叠后地震数据体。图3(c)是应用本发明的反褶积叠后地震数据体，其中初始给定的空间正则化参数均为0.1，井的约束为真实模型0-5Hz的信息。由图可知，无正则化项的反褶积叠后地震数据体是明显不稳定的；而本发明的反褶积叠后地震数据体是稳定的，且低频信息被较好地恢复。此外，相比于原始资料，应用本发明得到的叠后地震数据体的纵向分辨率得到了改善，且横向信息具有更好的连续性。
[0073] 图4示出了在另一应用实例中应用本发明得到的A井、B井、C井的井控多道反褶积结果与井滤波反射系数对比图。
[0074] D气田地表为沙漠和半沙漠环境，地震资料的频率和信噪比均偏低，地震资料的分辨率也很低，目前三维地震资料的主频为25Hz～30Hz，有效高频也只是50Hz，10米左右或小于10米砂体的识别非常困难，因此提高地震资料分辨率处理必不可少。采用本发明具体过程如下：
[0075] (1)选定初始空间正则化参数为20，将多道原始地震记录进行空间反褶积，得到地震反射系数系列；
[0076] (2)对井曲线进行滤波处理后，提取井的反射系数系列；
[0077] (3)将井的反射系数系列与地震反射系数系列进行运算得到的误差作为约束；
[0078] (4)此时井的反射系数系列与地震反射系数系列误差较大；于是，修改空间正则化参数为25，重复上述过程；
[0079] (5)经迭代后，当误差小于设定阈值时(此时空间正则化参数为30)，输出反射系数系列，并将之与原始地震记录提取的稳定子波进行褶积并叠加，得到最终的高分辨率地震数据体。
[0080] 图5(a)是某工区内一条原始地震剖面，图5(b)是应用本发明得到的地震剖面(最终确定的空间正则化参数均为30)。可以发现井控多道反褶积方法保持了地震响应的横向稳定性、连续性，且有效提高了地震资料的分辨率。处理后的地震剖面反射特征丰富，符合地质认识的砂体相变快的特征，与实钻测井滤波后波阻抗对比，反褶积后与实钻结果吻合，且可以识别最薄储层为8米厚的储层。
[0081] 本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。
[0082] 以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易 见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
