本发明涉及石油勘探技术，特别是关于石油勘探的数据处理技术，具体的讲是一种多次波识别方法及装置。
背景技术：
在陆上地震勘探中，地震波通过地下岩层传播，当地下有如基岩面、不整合面、火成岩面等强反射界面时，会产生层间多次波。这类多次波与深层一次波振幅相当，动校正时差较小，与深层一次波相干涉，难以识别，影响地震数据偏移质量，给后续的解释工作带来困扰。目前层间多次波压制方法以叠前为主，多次覆盖技术以来，人们认识到地震数据中的多次波可以根据其时距曲线特征加以识别。一般情况下，地震波传播速度随着深度的增加而增加，与多次波相比，一次波传播速度要高。共中心点(cmp)道集中，多次波与一次波同相轴之间出现倾角或时间差异，nmo之后，一次波是平的，而多次波由于nmo校正时间不足是下拉的，通常可以采用一个数学变换的方法将地震数据变换到新的数据域，从而将一次波与多次波投影到可分离的不同区域。但是，如果地层速度发生反转，与浅层高速体有关的多次波速度会接近该高速体以下的一次波速度，使得两者干涉，难以分离。在很多地区，在叠前阶段经过多次波压制后，仍然会有较强能量、频带宽的多次波残留。在叠后阶段，可以通过叠加地震剖面上强反射界面和强反射界面之下地震反射的形态和时间识别多次波。当地下反射界面倾角较小时，如果在强反射界面下方存在另一组与其形态一致，且出现时间小于其2倍的反射同相轴时，就可以判断该组同相轴为多次反射波。对于倾斜地层，叠加剖面上二次反射波的倾角约为其一次反射波的2倍。现有技术中，叠后阶段识别多次波的方法相对单一，主要依靠地震剖面上地震反射的外观、形态判断，主观性较强，识别的结果说服力不足。当强反射界面较多、多次波来源可能性较多，多次波波场特征复杂时，无法难以判断其来源层位。技术实现要素：针对当前多次波识别方法的不足，本发明实施例提供一种多次波识别方法，包括：获取目标区域的叠后地震数据和测井数据；根据所述测井数据确定各反射界面的一次反射系数；将所述一次反射系数中大于预设的反射系数阈值的反射界面作为强反射界面；根据确定的强反射界面和一次反射系数确定多次反射系数；根据所述一次反射系数、多次反射系数及叠后地震数据确定合成地震记录；所述合成地震记录包括：一次波的正演记录和包含多次波的正演记录；在所述合成地震记录中创建时窗窗口，根据各时窗中实际地震道与合成地震记录中各正演记录的相关系数识别多次波。本发明实施例中，所述的测井数据包括：声波时差曲线、密度测井曲线以及深度数据。本发明实施例中，所述的根据所述测井数据确定目标区域的一次反射系数包括：根据声波时差曲线确定速度数据；根据速度数据和深度数据确定双程旅行时间；根据速度和密度曲线确定波阻抗；对所述波阻抗按等时时间间隔dt方波化处理确定方波化的波阻抗；根据确定的方波化阻抗和下式确定一次反射系数；其中，rj为一次反射系数；z1为方波化的波阻抗；j为采样点序号，z1j+1为第j+1个采样点的方波化阻抗，z1j第j个采样点的方波化阻抗。本发明实施例中，所述的根据确定的强反射界面和一次反射系数确定多次反射系数包括：以等时时间间隔dt划分地层；利用反射率法确定划分的各地层的上行波能量和下行波能量，将各地层的上行波能量和下行波能量的比值作为全井段包含多次波的反射系数；利用所述强反射界面筛选所述全井段包含多次波的反射系数确定多次反射系数。本发明实施例中，所述的根据所述一次反射系数、多次反射系数及叠后地震数据确定合成地震记录包括：从所述叠后地震数据中提取统计子波；统计子波与反射系数序列褶积确定正演记录；其中，所述反射系数序列为：一次反射系数、全井段包含多次波的反射系数以及多次反射系数组成的反射系数序列；将所述正演记录进行归一化处理生成合成地震记录，所述合成地震记录包括：一次波正演记录sy，包含多次波的正演记录syn，其中n＝0，1...n。本发明实施例中，所述的在所述合成地震记录中创建时窗窗口，根据各时窗中实际地震道与合成地震记录中各正演记录的相关系数识别多次波包括：在合成地震记录中创建时间窗口；确定各时间窗口中，实际地震道的正演记录分别与一次波正演记录sy、包含多次波的正演记录sy0的相关系数；实际地震道的正演记录与正演记录sy0的相关系数大于实际地震道的正演记录与正演记录sy时，根据多次波的正演记录syn中波组特征的变化确定产生多次波的地层。本发明提供一种剥层法分析和识别多次波的方法和装置，利用反射率法正演模拟，了解地震波传播特点、识别多次波波场信息。通过此方法可以直接针对某些复杂层位，识别出产生的多次波的地层，在常规的叠加剖面上进一步识别出残留的多次波。为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1为本发明公开的多次波识别方法的流程图；图2为本发明实施例中的流程图；图3a为过a井0～4s的叠后地震剖面；图3b为本发明实施例中a井0～5800m的声波时差曲线ac；图3c为本发明实施例中a井0～5800m的密度曲线density；图4为本发明实施例中地震道及地震记录曲线；图5为本发明实施例中rj层的地震波反射路径；图6为本发明实施例中rj层的地震波反射路径；图7本发明实施例中syn中1.2s到2.4s的地震记录。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。本发明提供一种多次波识别方法，如图1所示，该方法包括：步骤s101，获取目标区域的叠后地震数据和测井数据；步骤s102，根据所述测井数据确定各反射界面的一次反射系数；步骤s103，将所述一次反射系数中大于预设的反射系数阈值的反射界面作为强反射界面；步骤s104，根据确定的强反射界面和一次反射系数确定多次反射系数；步骤s105，根据所述一次反射系数、多次反射系数及叠后地震数据确定合成地震记录；所述合成地震记录包括：一次波的正演记录和包含多次波的正演记录；步骤s106，在所述合成地震记录中创建时窗窗口，根据各时窗中实际地震道与合成地震记录中各正演记录的相关系数识别多次波。本发明在于通过剥层法识别和分析多次波，通过正演剥层前后含有多次波的正演记录，与一次波的正演记录、实际地震道对比分析，识别该层可能产生的多次波，最后识别出实际地震数据中的残留多次波和多次波来源。本发明实施例识别多次波的步骤，主要包括：获取地震数据、声波时差、密度测井曲线。根据声波时差求取倒数得到速度，并对测井曲线进行深时转换；速度和密度相乘得到波阻抗；对波阻抗激进行等时间隔方波化处理；由反射系数公式计算得到一次反射系数序列。根据一次反射系数序列的大小排序筛选出强反射地层。利用获得的一次反射系数序列通过反射率法计算得到全井段包含多次波的反射系数序列；结合筛选的强反射界面，依次去除强反射层及其以上的反射系数，重新通过反射率法计算得到强反射层以下井段包含多次波的反射系数序列。一次反射系数序列、包含多次波的反射系数序列与子波褶积得到正演记录；从叠后剖面提取测井位置处的实际地震道；对实际地震道、正演记录振幅归一化处理。创建时窗，分析时窗内实际地震道和正演记录之间的波组特征和互相关系数识别出实际叠后剖面中的残留多次波和来源。图2为本发明实施例的流程图，如图2所示，该方法包括：s201，获取地震数据和测井数据，具体为：1)获取叠后地震数据2)获取测井资料，包括声波时差和密度曲线如图3a为过a井0～4s的叠后地震剖面，图3b、图3c分别为a井0～5800m的声波时差曲线ac和密度曲线density。s202，计算一次反射系数：1)测井曲线深时转换：声波时差曲线ac求取倒数得到速度vp，并利用深度depth和速度vp关系求得双程旅行时间。过程如下：其中，n记录深度depth数据个数，time记录双程旅行时间数据。2)计算方波化波阻抗曲线：由于声波和密度测井测量方法和仪器所限，测井曲线会包含较多异常值，存在较多高频噪声。因此需要对测井曲线进行方波化处理，消除薄层和异常值的干扰，保持厚层界面信息，使测井曲线频带接近地震频带。速度vp和密度density相乘计算得到波阻抗。z＝vp×density对波阻抗z按等时时间间隔dt进行方波化处理得到的波阻抗z1，并得到方波化的时间time1。其中n为方波化时间间隔dt对应的方波化长度，z1为方波化波阻抗；dt、ρ分别为地层声波时差和密度；n为方波化间隔dt对应的方波化长度。如图4所示，a为波阻抗z，b为方波化的波阻抗z1。本实例中的方波化间隔dt为0.004s，过程如下：其中，m记录方波化个数，mean为方波化函数。3)计算一次反射系数：图4中的c为由公式计算得到的一次反射系数r，z1为方波化的波阻抗；j为采样点序号。程序实现过程如下：fori＝1:m-2r(i，1)＝(z1(i+1)-z1(i))/(z1(i+1)+z1(i))；ends203，筛选强反射界面：地层一次反射系数r从大到小排序并记录对应方波化的时间值time1。本实例中的排序利用matlab中的sort函数。[rrindex]＝sort(r(:，1))；time1(index)；其中rr为排序后的反射系数，time1(index)为对应的时间。如表1所示，本实施例中a井地层反射系数从大到小依次为0.2432、0.2211、0.2024、0.1918、0.1337，……，具体的，对应的强反射源界面time1时间值分别为1.9160s对应奥陶系底o1、1.7800s对应上二叠底p2l、1.1240s对应上三叠系须家河底t3x1、0.8000s对应侏罗系底j1、1.5040s对应下三叠系飞四段底t1f4、2.1400s对应寒武系筇竹寺底e1q。表1序号rtime1(s)10.24321.916020.22111.780030.20241.124040.19180.800050.16621.504060.13372.1400………………本发明实施例中，筛选出大于2/3最大反射系数的界面作为强反射界面，从而得到了5个强反射界面，按time1从小到大排序分别是：t1＝0.8000s，t2＝1.1240s，t3＝1.5040s，t4＝1.7800s，t5＝1.9160s。s204，计算多次反射系数：以方波化时间间隔dt为间隔，将地层划分为m层(方波化个数)，每层用time1表示。图5、图6为rj层(对应时间time1(j)层)的地震波反射路径，rj表示反射界面的一次波下行反射系数，利用反射率法求多次波的反射系数即求取图中rj层的地表的上行波能量u。一般地，对于任意反射系数为rj的反射层，下行波能量的公式为：dj+1＝-rj*u+(1-rj)*dj上行波能量的公式为：u＝rj*dj+(1+rj)*u由此得到每个time1层的上行波能量u，将上行波能量与下行波能量的比值作为多次反射系数。如图3的d到i所示，根据一次反射系数r由反射率法计算得到全井段包含多次波的反射系数r0；结合筛选强反射界面t1，去除t1层及其以上的一次反射系数r，重新由反射率法计算得到t1层以下井段包含多次波的反射系数r1，同理依次去除t2、t3、t4和t5层及其以上的反射系数分别得到r2，r3，r4，r5。具体实现过程如下：测井数据：深度声波时差密度都是数组；声波时差dt求倒数得到速度速度vp乘以密度ρ得到：波阻抗速度vp和深度depth得到时间每个时间点i都对应一个波阻抗i＝1～n；方波化：因为不是等时间隔采样，time1表示的是每层的包括所有层，以dt为时间间隔重采样得到：每个时间点j都对应一个波阻抗j＝1～m；由公式，通过z1计算一次反射系数得到进一步，关于多次反射系数计算：更新上行波公式：u＝rj*dj+(1+rj)*u；更新下行波公式：dj+1＝-rj*u+(1-rj)*dj；初始条件：下行波能量为1，上行波能量为0，即d1＝1,u＝0；更新上行波公式：u＝r1*d1+(1+r1)*u＝r1依次计算下去得到整个多次反射系数r1的计算：比如本实施例中，一次反射系数中r010为最大值，是最大的反射界面，去掉r010及以上的数据得到：这时只有m-10个数据了。多次反射系数初始条件：下行波能量为1，上行波能量为0，即d1＝1,u＝0；更新上行波公式：u＝r11*d1+(1+r11)*u＝r11依次计算下去得到整个代码过程如下：其中，startt为起始计算时间。s205，计算合成地震记录：1)提取地震子波，合成地震记录：以方波化时间间隔dt为采样间隔，从叠后地震数据中提取a井处实际地震道s，并从地震道s中提取出统计子波w；反射系数序列r、r0、r1、r2、r3、r4、r5与地震子波w褶积得到的正演记录sy、sy0、sy1、sy2、sy3、sy4、sy5。syn＝rn*w(n＝，0，1，2，3，4，5)syn为测井合成地震记录；w为地震子波；rn为反射系数。2)地震道能量归一化处理：将实际地震道s及正演记录sy、sy0、sy1、sy2、sy3、sy4、sy5的振幅分别除以该道振幅绝对值的最大值，更新原记录。syn＝[s，syn]syn＝syn./max(abs(syn))本实例中，提取实际地震道、统计子波借用了地震反演软件hrs的strata模块，褶积运算借用了matlab中的conv函数，过程如下：loads.txt；％加载实际地震道syn＝[s，syn]；syn＝syn./max(abs(syn))；如图7所示，为syn中1.2s到2.4s的地震记录；s106，识别多次波：1)确定多次波时窗：在地震记录syn中，创建时间窗口，对比时间窗口中一次波的正演记录sy和实际地震道s的波组特征，并求互相关系数。理论上，一次波的正演记录sy和实际地震道s二者波组特征越相近，互相关系数越大，地震道多次波干扰越少。反之，二者差异越大，互相关系数越小，多次波干扰越严重，甚至一次波能量被多次波掩盖。因此，二者匹配不好相关性差的地方说明地震数据中多次波发育。如图7所示，1.2s～2.0s的时间范围内，一次波合成地震记录sy与实际地震道s符合很好，说明该段实际地震道以一次有效波为主。2.0s～2.4s的时间范围内，一次波合成地震记录sy与实际地震道s波组特征有较大差异，而含有多次波的正演记录sy0与实际地震道s符合更好。具体表现在强反射振幅的个数不同，实际地震道s、含有多次波的正演记录sy0明显有多个强能量振幅。说明2.0s～2.4s内，除一次波反射之外，还包含有较强的多个多次波反射，实际地震资料在该段时间范围内可能残留了较强的多次波干扰。在本实例1.2s到2.4s的地震记录中，以0.2s为间隔创建6个时间窗口，并求每个窗口中实际地震道s与一次波合成地震记录sy、包含多次波的正演记录sy0的互相关系数。本实例中的相关运算借用了matlab中的corrcoef函数。如表2所示，2.0～2.2s，2.2～2.4s的时间窗口，一次波正演记录sy和实际地震道s的互相关系数明显偏小，反观包含多次波的正演记录sy0和实际地震道s的互相关系数更大，说明在2.0～2.2s，2.2～2.4s的时间窗口内，多次波干扰越严重。表2综上，在本实施例中，实际地震数据中多次波时窗范围为：2.0～2.2s，2.2～2.4s。2)确定产生多次波的来源地层：识别多次波特征在于，在同一时窗内，对比合成地震记录和实际地震道的波组特征和相关性。理论上，剥层前包含多次波的正演记录和实际地震道二者波组特征更相近，而剥层后包含多次波的正演记录和一次合成地震道二者波组特征更相近。理论上，在同一时窗内，剥层前包含多次波的正演记录和实际地震道s二者波组特征更相近，而剥层后包含多次波的正演记录和一次合成地震道sy0二者波组特征更相近。如图4所示，图4的a为采用本发明实施例的a井处实际地震道s，图4的b到h为采用本发明实施例的a井的合成地震记录。在2.0～2.2s时窗范围内2.1s附近，含有多次波的合成地震记录sy1较sy0，波组特征上看无明显变化，说明该层产生的多次波相对较弱，不是主要的多次波来源层，而sy2较sy1变化却比较明显，到sy3、sy4、sy5已无明显变化，且接近一次波合成地震记录s的波组特征。说明t2层位是产生2.1s处多次波的主要地层。同理，在2.2～2.4s时窗范围内2.3s附近的多次波来源于t3层，2.35s附近的多次波来源于t1层。依次求取每个时窗中，正演记录sy0、sy1、sy2、sy3、sy4、sy5与实际地震道s、一次合成地震道sy的互相关系数。如表3所示，随着剥层的进行，正演记录sy0、sy1、sy2、sy3、sy4、sy5与实际地震道s的互相关系数在减小，而与一次合成地震道sy的互相关系数在增大，剥去t2层后corrcoef(s，sy2)明显减小，而corrcoef(sy，sy2)明显增大，说明2.0～2.2s窗口中的多次波来源于t2层。同理，如表4所示，2.2～2.4s窗口中的多次波主要来源于t1，t3层。表32.0～2.2s时窗内的互相系数表表42.2～2.4s时窗内的互相系数表corrcoef(，)ssysy00.84200.7178sy10.75320.7832sy20.76020.8375sy30.70360.9110sy40.68090.9544sy50.64220.9886综上，在本实施例中，2.0～2.2s时窗内的多次波主要来源于t2地层，即上三叠系须家河底t3x1。2.2～2.4s时窗内的多次波主要来源于t1和t3地层，即侏罗系底j1和下三叠系飞四段底t1f4。本发明可适用于地震解释工作中，识别中深层具有强反射界面的地区产生的层间多次波。这些强反射界面通常是基岩面、不整合面、火成岩面等，反射波在返回地面遇到这些界面时，可能再次反射后向下传播，形成多次波，这些多次波与深层一次波能量相当，混在一起难以分离。随着勘探的深入，更深深度地层逐渐进入开发阶段，由于深层多层波的干扰，面临高产井地震响应不明确、确立高产井模式困难、有效储层预测的问题。本发明可分析产生层间多次波的来源层位，使得多次波的压制针对性更强，一定程度上减少了压制多次波的不确定性，更易分离和消除多次波干扰。本发明通过反射率法正演模拟和波组对比，剥层分析和识别多次波。目前叠后阶段识别多次波的方法相对单一，主要依靠地震剖面上地震反射的外观、形态判断，主观性较强，多解性强，识别的结果说服力不足。本发明采用和测井、地震解释相结合的方法，能实现半定量的多次波识别，从而大大减少了多次波识别的不确定性，更好地分离和消除多次波干扰波。本发明实施例中，一次波和多次波合成地震记录正演，通过计算合成地震记录和实际地震道的相关系数，实现半定量的多次波识别；本发明可以直接针对某些复杂层位，识别出产生多次波的地层，极大地减少无关数据的分析和运算。应当注意，尽管在附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。虽然本申请提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的手段可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的装置或客户端产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行(例如并行处理器或者多线程处理的环境，甚至为分布式数据处理环境)。术语“包括””、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、产品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、产品或者设备所固有的要素。在没有更多限制的情况下，并不排除在包括所述要素的过程、方法、产品或者设备中还存在另外的相同或等同要素。上述实施例阐明的单元、装置或模块等，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申请时可以把各模块的功能在同一个或多个软件和/或硬件中实现，也可以将实现同一功能的模块由多个子模块或子单元的组合实现等。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。本领域技术人员也知道，除了以纯计算机可读程序代码方式实现控制器以外，完全可以通过将方法步骤进行逻辑编程来使得控制器以逻辑门、开关、专用集成电路、可编程逻辑控制器和嵌入微控制器等的形式来实现相同功能。因此这种控制器可以被认为是一种硬件部件，而对其内部包括的用于实现各种功能的装置也可以视为硬件部件内的结构。或者甚至，可以将用于实现各种功能的装置视为既可以是实现方法的软件模块又可以是硬件部件内的结构。本申请可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构、类等等。也可以在分布式计算环境中实践本申请，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本申请可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如rom/ram、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，移动终端，服务器，或者网络设备等)执行本申请各个实施例或者实施例的某些部分所述的方法。本说明书中的各个实施例采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。本申请可用于众多通用或专用的计算机系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的电子设备、网络pc、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。当前第1页12