CN112946742A - 一种拾取精确叠加速度谱的方法 - Google Patents

一种拾取精确叠加速度谱的方法 Download PDF

Info

Publication number
CN112946742A
CN112946742A CN202110284488.5A CN202110284488A CN112946742A CN 112946742 A CN112946742 A CN 112946742A CN 202110284488 A CN202110284488 A CN 202110284488A CN 112946742 A CN112946742 A CN 112946742A
Authority
CN
China
Prior art keywords
data
velocity
gather
spectrum
frequency
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.)
Granted
Application number
CN202110284488.5A
Other languages
English (en)
Other versions
CN112946742B (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.)
Chengdu Jiekesi Petroleum Natural Gas Technology Development Co ltd
Sichuan Changning Natural Gas Development Co ltd
Original Assignee
Chengdu Jiekesi Petroleum Natural Gas Technology Development Co ltd
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 Chengdu Jiekesi Petroleum Natural Gas Technology Development Co ltd filed Critical Chengdu Jiekesi Petroleum Natural Gas Technology Development Co ltd
Priority to CN202110284488.5A priority Critical patent/CN112946742B/zh
Publication of CN112946742A publication Critical patent/CN112946742A/zh
Application granted granted Critical
Publication of CN112946742B publication Critical patent/CN112946742B/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. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • 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. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • 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. for interpretation or for event detection
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种拾取精确叠加速度谱的方法，包括以下步骤S1,确定相关的叠加道上层位数据,根据设计的分频参数对道集数据进行滤波处理，提取速度分析点上的分频道集数据体；S2，对速度分析点中不同的分频道集数据体实施各个层位数据的引入，得到相关速度分析点上的各个分频道集上的层位数据；S3，以相关速度分析点上的各个分频道集上的层位数据、初始速度谱作为约束，依次对各个速度分析点上的各个分频道集数据体进行相关速度谱的迭代性求取，得到各个速度分析点上的叠加速度谱；本发明解决了现有技术存在的拾取的速度谱精度不高等问题。

Description

一种拾取精确叠加速度谱的方法
技术领域
本发明涉及石油天然气勘探信号处理技术领域，具体是一种拾取精确叠加速度谱的方法。
背景技术
随着页岩气勘探技术的发展，近年来，四川盆地内海相页岩气勘探及开发取得巨大进展。如四川盆地的长宁-威远区块、威荣页岩气勘探区及焦石坝、平桥页岩气勘探区块中，大批页岩气水平井经压裂后获得工业气流，取得可观的经济效益。在对海相页岩气勘探及开发实践中，大量勘探成果表明良好的保存条件是页岩气高产富集的关键因素之一。而在靠近大断裂带附近的页岩气勘探区，由于构造复杂，造成地震反射成像困难，影响到相关的页岩气井布设及勘探进展。
在地震资料处理过程中，速度分析是关键技术之一。速度分析的技术流程是以不同的速度对共中心点道集进行扫描，以动校正叠加后能量团的聚焦性作为拾取标准，因此高精度的动校正技术是影响速度谱精度的主要因素之一。动校正的目的是消除炮检距对反射波旅行时的影响，将共中心点道集(CMP道集)中的反射波同相轴拉平，增强利用叠加技术压制干扰的能力，减小叠加过程引起的反射波同相轴畸变，从而提高地震数据的信噪比。其精确性直接影响干扰波的压制效果以及后续偏移成像的质量。因此，速度分析在地震资料成像中是相当关键的一个步骤。
由于速度谱的拾取对地震资料的成像及后续处理相当关键，因此，相关地球物理工作者对其拾取的精确性相当关注。一些发明专利如CN201810571640.6的发明专利《一种基于深度强化学习的叠加速度谱拾取方法及处理终端》主要步骤为获取包括地震反射波的原始共中心点地震道集数据并计算出以各个时刻下的最优扫描速度组成的叠加速度谱，将叠加速度谱输入至自编码网络，得到编码后的高阶能量团特征，将高阶能量团特征编码输入至策略网络，拾取各个时刻下的最优扫描速度，输出速度序列，对速度序列进行评价，输出奖励值，根据奖励值来训练策略网络；如专利号CN201610581659.X《基于速度模型的速度谱解释方法》利用地面地震叠前数据共中心点(CMP)道集数据产生速度谱，确定初始速度模型函数，基于初始速度模型函数，在CMP中选择至少一个种子点，求取种子点处的初始速度模型函数，以种子点处的初始速度模型函数作为约束，在速度谱上沿时间轴拾取谱能量最大的点的速度。但是，针对速度谱精确拾取和分析的技术还不够完善，现有技术主要是利用道集上一个设定速度(从上一个速度谱分析拷贝出来)，对道集进行基于拷贝速度进行层拉平后的速度分析，道集由于没有分频处理，一些高频响应的能量团受到压制，不利于精确的速度分析；而且没有层位数据进行约束，速度分析的道集上没有层位信息。
主要体现在以下方面：
(1)采用常规道集资料进行速度谱的拾取，由于速度分析所使用的道集中存在不同的频段，往往由于遭受干扰等导致拾取的速度谱精度不高。
(2)速度分析道集中不同时间范围内，由于大地滤波作用及叠加次数等因素，导致道集的上、中、下区域的信噪比及一些目的层响应受到压制，也存在一些噪音干扰。
发明内容
为克服现有技术的不足，本发明提供了一种拾取精确叠加速度谱的方法，解决现有技术存在的拾取的速度谱精度不高等问题。
本发明解决上述问题所采用的技术方案是：
一种拾取精确叠加速度谱的方法，包括以下步骤：
S1,确定相关的叠加道上层位数据,根据分频参数对道集数据进行滤波处理，提取速度分析点上的分频道集数据体；
S2，对速度分析点中不同的分频道集数据体实施各个层位数据的引入，得到相关速度分析点上的各个分频道集上的层位数据；
S3，以相关速度分析点上的各个分频道集上的层位数据、初始速度谱作为约束，依次对各个速度分析点上各个分频道集数据体进行迭代性的速度谱的求取，得到各个速度分析点上的叠加速度谱。
作为一种优选的技术方案，还包括位于步骤S3之后的以下步骤：
S4，判断各个速度分析网格点的叠加速度数据的收敛情况，若未收敛，则返回执行步骤S1～S3；若收敛，则结束处理。
作为一种优选的技术方案，步骤S1包括以下步骤：
S11，准备地质、测井以及地震数据；
S12，建立速度谱拾取约束使用的各个层位数据；
S13，设定关于叠前道集的分频处理参数。
作为一种优选的技术方案，步骤S11中，采用带通滤波器对道集数据进行滤波处理。
作为一种优选的技术方案，步骤S2包括以下步骤：
S21,基于叠后处理资料中所确定的初始速度谱及速度分析点叠加道上各个层位数据，确定其相关层位的地震反射特征及双程反射时间数据；
S22,对相关速度分析点上的整个频段道集的相关层位数据进行引入，得到该速度分析点道集上的层位数据；
S23，根据速度分析点上的道集上各个层位数据信息提取其分频道集数据，利用道集上的各个层位数据信息对其分频道集数据进行相关层位数据投影，得到相关分频道集的层位数据。
作为一种优选的技术方案，所述的层位数据信息包括但不限于层位的双程反射时间、坐标信息、偏移距信息。
作为一种优选的技术方案，步骤S3包括以下步骤：
S31,对各个速度分析点上的各个分频道集数据利用相关的层位数据、初始速度谱作为约束，对相关的分频体实施有顺序性、迭代性的速度谱拾取，得到各个速度分析点上的叠加速度谱；
S32,利用步骤S31获得的各个速度分析点上的叠加速度谱进行叠加处理，得到相关的叠加数据体，再对相关目的层位进行剩余静校正量的计算，并利用剩余静校正量对道集进行校正。
作为一种优选的技术方案，步骤S31中对相关的分频体实施从低频段向高频段的排列顺序的速度谱拾取。
作为一种优选的技术方案，步骤S31中还包括以下步骤：对中低频段、低频段所获得的叠加速度谱在双程反射时间上进行垂向上的速度点之间的加密。
本发明相比于现有技术，具有以下有益效果：
(1)本发明设定速度分析网格点，并对相关分频道集进行关于网格点的道集重构后，指导对不同分频道集数据的各速度分析网格点的速度拾取，按频段的频率从低到高拾取速度谱并进行叠加，这样速度拾取的精度越来越高；避免了由于遭受干扰等导致拾取的速度谱精度不高的问题；本发明可以实现精确的叠加速度分析，从而提高地震资料的成像精度，进而减少钻井风险，提高相关油气勘探的经济效益；
(2)本发明设计一系列层位数据，并对道集进行投影，提示对这些层位所对应的道集上反射层进行精确的速度拾取后实现层拉平，有层位数据进行约束，速度分析的道集上有层位信息；并且，对相关层位上所对应的各个分频道集上能量团聚焦的情况进行分析，也能判断分频参数的合理性；通过实现精确的叠加速度分析，并能经过相关计算获得准确的层速度数据体，从而提高地震资料的成像精度及储层预测的准确性，进而减少钻井风险，提高相关油气勘探的经济效益；
(3)本发明频段分布参数的设定，可以针对速度分析网格点上的整个道集进行，也可以相关目的层位分别设定各自的分频参数进行处理；
(4)本发明设定解释网格后对相关界面进行追踪解释，从而得到相关反射层界面的层位数据；对叠后偏移地震数据体的层位进行内插后，再根据相关层位上的层位数据信息(如双程反射时间、坐标信息)对三维叠加数据体进行投影，得到相关叠加数据体的层位数据，再对其进行相关速度分析网格点上叠加道上的层位数据的提取；设定关于叠前道集的频段分布参数，对道集设定一系列的滤波参数，对三维地震道集数据体分别进行所设定的一系列滤波参数处理，得到一系列滤波范围从低到高的道集数据体；
(5)本发明频段分布参数的设计主要是在处理频段内对相关几个频率滤波范围的设计，从而得到相关的滤波数据，相关频率滤波范围可以设定为等间隔或不等间隔；
(6)本发明对速度分析点上的整个频段道集(没有进行分频处理的道集)上的各个层位数据的获得，主要操作为根据相关速度分析点上叠加道的层位数据的双程反射时间和地震反射特征，以某一层位为例，将其以某一入射角度的道集数据上的某一采样点位置的投影，原则上要求投影点的地震反射特征与双程反射时间与解释的层位数据的双程反射时间和地震反射特征相吻合；并根据动校正处理的道集(初始速度谱处理)，根据确定的固定时窗在该道集内进行该层位的层位自动追踪，从而得到该层的层位数据。依次类推，完成各个层位数据在该速度分析点上的道集内层位数据的确定。并对叠加道与各个入射角地震道进行关于目的层的相关系数计算，得到一系列的相关系数数据值；从而得到相关的某一入射角度的道集数据上的层位数据的投影；
(7)本发明鉴于高频段的分频数据体中存在分辨率相对较高的信息，在对高频段的分频道集的速度谱拾取过程中，对低频和中频段所获得的叠加速度谱在双程反射时间上进行速度分析网格点垂向上之间的加密，从而使其相关的反射同相轴达到同相叠加的效果。并根据相关的叠加速度谱的迭代处理，获得更为准确的叠加速度谱；
(8)本发明在速度谱的拾取中，利用相关的层位数据对其地震反射波进行校平后，所拾取的速度谱点作为相关的速度谱数据点。依次类推，完成各个速度分析点上的精确速度谱点的拾取；各个速度分析点上的叠加速度谱进行叠加处理，得到相关的叠加数据体，再对相关目的层位进行剩余静校正量的计算，并利用剩余静校正量对道集进行校正后，再次进行相关速度分析点上的速度谱的精确拾取，并根据相关叠加速度数据收敛情况终止叠加速度的拾取，从而得到更为精确叠加速度谱。
附图说明
图1为本发明所述的一种拾取精确叠加速度谱的方法的流程图。
具体实施方式
下面结合实施例及附图，对本发明作进一步的详细说明，但本发明的实施方式不限于此。
实施例1
如图1所示，一种拾取精确叠加速度谱的方法，包括以下步骤：
S1,确定相关的叠加道上层位数据,根据分频参数对道集数据进行滤波处理，提取速度分析点上的分频道集数据体；
S2，对速度分析点中不同的分频道集数据体实施各个层位数据的引入，得到相关速度分析点上的各个分频道集上的层位数据；
S3，以相关速度分析点上的各个分频道集上的层位数据、初始速度谱作为约束，依次对各个速度分析点上各个分频道集数据体进行速度谱的求取，得到各个速度分析点上的叠加速度谱。
本发明设定速度分析网格点，并对相关分频道集进行关于网格点的道集重构后，指导对不同分频道集数据的各速度分析网格点的速度拾取，按频段的频率拾取速度谱并进行叠加，这样速度拾取的精度越来越高。避免了由于遭受干扰等导致拾取的速度谱精度不高的问题。值得说明的是，本发明的初始速度谱，可以从原始速度谱分析中得到。
本发明可以实现精确的叠加速度分析，从而提高地震资料的成像精度，进而减少钻井风险，提高相关油气勘探的经济效益。
利用本发明技术所得到的准确的速度谱成果，与采用常规速度谱分析结果所得到的地层层速度模型(经均方根速度计算所得到)进行相关后续的叠后时间偏移处理后，以两者过井的地震剖面进行对比分析。相关对比成果表明，本发明技术所得到的地震剖面优于常规处理所得到的地震剖面，这也证明了本发明技术对精确速度谱的拾取是有效的。
具体地，包括以下执行步骤：
(1)确定叠后偏移数据体的各个层位数据，经内插后得到各个速度分析点上的叠加道的层位数据，确定叠前动校正处理后的道集数据体关分频处理参数并经相关带通滤波处理后得到一系列的分频道集数据体；
(2)基于道集上的初始速度谱及相关速度分析点上的叠加道上各个层位数据并对相关分频道集数据进行投影，从而得到相关速度分析点上的各个分频道集上的层位数据。
(3)以相关速度分析点上的各个分频道集上的层位数据、初始速度谱进行约束，依次对各个速度分析点上各个分频道集数据体进行精确的迭代性速度谱的求取，从而得到各个速度分析点上的精确叠加速度谱。
作为一种优选的技术方案，还包括位于步骤S3之后的以下步骤：
S4，判断各个速度分析网格点上相关层位上叠加速度数据的收敛情况，若未收敛，则返回执行步骤S1～S3；若收敛，则结束处理。
利用迭代，可实现更加精确地速度谱拾取，更加便于分析，速度谱拾取和分析地准确率更高。
作为一种优选的技术方案，步骤S1包括以下步骤：
S11，准备地质、测井以及地震数据；
S12，建立速度谱拾取约束使用的各个层位数据、初始速度谱；
S13，设定关于叠前道集的分频处理参数。
作为一种优选的技术方案，步骤S11中，采用带通滤波器对道集数据进行滤波处理。
作为对步骤S1的详细实施和描述,步骤S1用于确定相关的叠加道上层位数据,根据设计的分频参数对道集数据进行滤波处理，提取速度分析点上的分频道集数据体。具体如下：
本发明的步骤S1中，确定叠后偏移数据体的各个层位数据，经内插后得到各个速度分析点上的叠加道的层位数据，确定叠前动校正处理后的道集数据体关分频处理参数并经相关带通滤波处理后得到一系列的分频道集数据体，具体包括以下步骤：
(1-1)地质、测井以及地震数据准备。地质数据包括岩芯录井资料、地质分层资料及岩石物理测试成果等，地震数据为常规三维叠前偏移道集、叠加数据体或叠后时间偏移地震数据体，叠后速度分析的初始速度谱拾取的相关CMP速度分析点及其对应的初始速度谱；
(1-2)建立速度谱拾取约束使用的各个层位数据。该层位数据是指利用叠后偏移地震数据体、测井资料、地质分层数据等实施相关井-震标定后，在井上确定相关有意义的反射层界面，设定解释网格后对相关界面进行追踪解释，从而得到相关反射层界面的层位数据。另外，对解释的层位数据进行内插、圆滑等处理后，得到1线X1道的层位数据。依次类推，完成各个层位的解释，并得到相关的层位数据。经上述操作后，得到各个速度分析点上叠加道上的各个层位数据。其中，相关速度分析点上的层位数据是指对各个层位进行内插、圆滑等处理后，得到相关速度分析点上的叠加道上的层位数据。其中，速度分析点的叠加道上层位数据的获取，其主要操作步骤为对叠后偏移地震数据体的层位进行内插后，再根据相关层位上的层位数据信息(如双程反射时间、x及y坐标)对三维叠加数据体进行投影，得到相关叠加数据体的层位数据，再对其进行相关速度分析点上叠加道上的层位数据的提取。
(1-3)设定关于叠前道集的分频处理参数。该步骤主要为对道集设定一系列的带通滤波参数，对三维地震道集数据体分别进行所设定的一系列带通滤波参数处理，得到一系列滤波范围从低到高的道集数据体。其中，带通滤波参数的设计具体操作为对目的层的地震主频进行分析及确定，并根据道集的整个处理频段范围，关于地震主频对处理的频段范围进行一系列的高、中、低频率段的设计，从而得到相关的分频处理参数。其中，带通滤波处理是公知技术，在本发明中不再详细阐述。如在实例中，应用离散傅里叶变换(DFT)把叠前道集数据变换到频率域，按设计的分频参数利用频率域滤波方法生成一系列分频数据体。此外，所述的叠前道集数据是指经过一系列叠前去噪、野外静校正处理及振幅补偿、反褶积、速度分析及剩余静校正迭代处理后的道集，也就是对最终叠加时所使用的道集数据，是利用叠前地震数据共中心点(CMP)道集。
优选地，地质、测井以及地震数据准备通过地球物理勘探、测井资料或者查表获取，这些数据主要是为了对井中相关层位进行标定及解释做准备。其中，初始速度谱是指经过n次常规速度分析及剩余静校正后的道集所拾取的速度谱，其所对应的道集(n次剩余静校正处理后，一般为3次)参与后续的分频处理计算。一般情况下，所谓的CMP速度分析点是指设定速度分析网格上的速度分析点，速度网格为n线Xm道(如20线X20道)，也可以在异常体上对其进行网格加密进行速度谱的拾取。原则上，层位是地震反射连续性相对较好，能量较强的同相轴所对应的地质层位。
优选地，道集的整个处理频段范围是指在得到最终叠加前的道集处理频率范围。所设计的一系列的高、中、低频率段，可以根据对目的层的成像效果、专家经验及计算精度等确定其范围及数目。一般情况下，所设定的滤波频段范围越小，则其数目相对较多，从而导致其速度分析点上的拾取工作量增大，但计算精度相对较高；反之，则计算精度相对较差，拾取工作量相对减小。其中，整个处理频段的范围可以从相关的处理流程参数中获得。原则上，要求经分频处理后的各个分频数据体在相关层位上的所对应速度分析的能量团能相对聚焦，或者能相对清楚地识别出速度谱上的层位上所对应的谱能量最大时所对应的速度值的位置；如在相关分频体上，某一个测试点上所对应的某一层位上的谱能量最大的点的速度值在从低频到高频上能逐渐收敛，相邻层位之间的一些能量团能在高频段中也相对于比低频段中能显示及识别出来。一般情况下，滤波频段范围及数目所用的测试点越多，则其更容易确定其分频参数，计算精度也越高；反之，会影响到分频参数的精确设计，导致其叠加速度的计算精度相对降低。在实际操作中，测试点的位置选择及数目主要是根据专家经验、速度分析精度及工作量等确定。一般情况下，测试点位置的选择应该选择在信噪比相对较高的区域，点与点之间的间隔应该大于或等于1000个CDP距离。
优选地，对叠前道集的分频处理，可以针对速度分析点上的整个道集进行，也可以相关目的层段分别设定各个速度分析点各自的分频参数进行处理，或者单独根据实际情况，对各个速度分析点设计其相关的分频处理参数。其中，目的层段是指两个层位之间的层段，可以包含一个及一个以上的层位。如何划分不同的目的层段，可以根据实际情况、专家经验及预测精度等确定。
优选地，由于道集的信噪比不高的时候，可以使用超级道集进行后续的计算。超级道集就是几个面元组成的道集，这个是公知技术。
优选地，分频处理参数的设计主要是在处理频段内对相关几个频率滤波范围的设计，从而得到相关的滤波数据。原则上，相关频率滤波范围可以设定为等间隔或不等间隔，也可以根据专家经验、实际地震资料及速度分析精度等进行设定。
一般情况下，包括起始频率
Figure BDA0002979869960000111
止频率
Figure BDA0002979869960000112
中心频率
Figure BDA0002979869960000113
频段长度Δdi、正整数i，相关的计算公式为：
Figure BDA0002979869960000114
Figure BDA0002979869960000115
其中，式中，i为每个频率区段的标号，
Figure BDA0002979869960000116
为滤波后第i频段起频率，
Figure BDA0002979869960000117
为滤波后第i频段止频率，
Figure BDA0002979869960000118
为第i频段的中心频率值，Δdi为第i频段长度。
当Δdi恒等于0，则为单频体。设计的
Figure BDA0002979869960000121
应有一些中心频率值与研究区内不同时间范围内(如从上至下)某一些目的层的调谐频率一致或相近，最小的起频率及最大的止频率应包括在有效频带范围内(整个处理频段)，频段与频段之间的可以重叠，但重叠范围应不大于0.25di。
优选地，所述Δdi值固定，Δdi＝△d，所述频段分布参数还包括
Figure BDA0002979869960000122
正整数k、正整数n，式中，k为频段数，
Figure BDA0002979869960000123
为设计的最后的第k个频段的中心频率，
Figure BDA0002979869960000124
为第一个频段的中心频率，n为步长，
Figure BDA0002979869960000125
此设计适用于设置对整个频段范围设计等间隔的频段分频情形。
优选地，举例对整个频段范围设计等间隔的分频参数，其主要操作为利用△d(设定为等间隔的固定值)和
Figure BDA0002979869960000126
两个参数，对整个频段范围设计其分频数目及步长，计算公式如下：
Figure BDA0002979869960000127
Figure BDA0002979869960000128
式中，k为频段数，
Figure BDA0002979869960000129
为设计的最后的第k个频段的中心频率，
Figure BDA00029798699600001210
为第一个频段的中心频率，n为步长，
Figure BDA00029798699600001211
此设计适用于设置对整个频段范围设计等间隔的频段分频情形。本发明设定速度分析网格点，并对相关分频道集进行关于网格点的道集重构后，再根据道集上层位信息，指导对不同分频道集数据的该网格点的速度拾取，按频段的频率从低到高拾取速度谱并进行叠加，这样速度拾取的精度越来越高，避免了由于遭受干扰等导致拾取的速度谱精度不高的问题。而且，本发明设计一系列层位数据，并对道集进行投影，提示对这些层位所对应的道集上反射层进行精确的速度拾取后实现层拉平，层位数量可以设计30多层，有层位数据进行约束，速度分析的道集上有层位信息。
本发明解决了速度分析道集中不同时间范围内，由于大地滤波作用及叠加次数等因素，导致道集的上、中、下区域的信噪比及一些目的层响应受到压制，以及存在一些噪音干扰导致拾取的速度谱精度不高的问题。
具体地，本发明可以首先针对密度计算相关的反演、属性数据体及建立模型用的层位数据；然后对不同层间建立相关的密度计算模型及优选属性组合，并经相关计算后得到不同层间的密度数据体；利用相关的波阻抗反演或属性计算方法实施计算后，得到波阻抗数据体，将波阻抗数据体及密度数据体代入相关的层速度计算模型进行计算后，得到一个层速度数据体，利用该数据体可以参与后续的层速度深度模型的建立，并实施叠前深度偏移处理。
以低频、中频、高频三个频段区分设置为例，先对低频道集上的层位信息，粗略地拾取相关层位上的所对应聚焦后的能量团，能使道集上的对应层位反射层进行拉平，从而得到低频道集上的速度谱；再利用这个速度谱对中频道集进行相关层位的反射层精确拉平，得到第二次速度分析的速度谱；再利用高频道集信息，对第二次速度谱再次精确分析，对高频道集上的相关层位的反射层进行精确拉平后，得到第三次速度谱，从而得到精确的速度谱。
本发明可以实现精确的叠加速度分析，并能经过相关计算获得准确的层速度数据体，从而提高地震资料的成像精度及储层预测的准确性，进而减少钻井风险，提高相关油气勘探的经济效益。
作为一种优选的技术方案，步骤S2包括以下步骤：
S21,基于叠后处理资料中所确定的初始速度谱及速度分析点叠加道上各个层位数据，确定其相关层位的地震反射特征及双程反射时间数据；
S22,对相关速度分析点上的整个频段道集的相关层位数据进行引入，得到该速度分析点道集上的层位数据；
S23，根据速度分析点上的道集上各个层位数据信息提取其分频道集数据，利用道集上的各个层位数据信息对其分频道集数据进行相关层位数据投影，得到相关分频道集的层位数据。
作为一种优选的技术方案，所述的层位数据信息包括但不限于层位的双程反射时间、坐标信息、偏移距信息。
作为对步骤S2的详细实施和描述,步骤S2用于对速度分析点中不同的分频道集数据体实施各个层位数据的引入，得到相关速度分析点上的各个分频道集上的层位数据。具体如下：
步骤S2基于道集上的初始速度谱及相关速度分析点上的叠加道上各个层位数据并对相关分频道集数据进行投影，从而得到相关速度分析点上的各个分频道集上的层位数据，包括以下步骤：
(2-1)基于叠后处理资料中所确定的初始速度谱及速度分析点叠加道上各个层位数据，确定其相关层位的地震反射特征及双程反射时间数据。在该步骤中，主要是根据对拾取初始速度谱中的速度分析点进行相关初始速度谱提取及叠加道上的层位数据信息。其中，相关的叠加道上的层位数据信息是指相关层位的地震反射特征及双程反射时间信息。
(2-2)对相关速度分析点上的整个频段道集的相关层位数据进行引入，得到该速度分析点道集上的层位数据。在该步骤中，具体操作为根据相关速度分析点上的初始速度谱作为背景数据，对该速度分析点上的整个频段道集数据进行针对初始速度谱的动校正处理；并根据该速度分析点上叠加道上各个层位数据的偏移距-双程反射时间数据对、坐标信息等对该点的道集数据相关数据道进行投影并展开自动追踪，从而得到该速度分析点上道集数据的各个层位数据。依此类推，完成各个速度分析点上道集的层位数据投影。
(2-3)根据速度分析点上的道集上各个层位数据信息并提取其分频道集数据，利用道集上的各个层位数据信息对其分频道集数据进行相关层位数据投影，得到相关分频道集的层位数据。其中，所述的层位数据信息是指层位的双程反射时间、x及y坐标信息、偏移距信息等。投影就是按相关的层位数据信息，将其投影到相关分频道集中相同信息的位置点上。
本发明的该步骤中，对速度分析点上的整个频段道集(没有进行分频处理的道集)上的各个层位数据的获得，主要操作为根据相关速度分析点上叠加道的层位数据的双程反射时间和地震反射特征，以某一层位为例，将其以某一入射角度的道集数据上的某一采样点位置的投影，原则上要求投影点的地震反射特征与双程反射时间与解释的层位数据的双程反射时间和地震反射特征相吻合；并根据动校正处理的道集(初始速度谱处理)，根据确定的固定时窗在该道集内进行该层位的层位自动追踪，从而得到该层的层位数据。依次类推，完成各个层位数据在该速度分析点上的道集内层位数据的确定。其中，所设定的固定时窗一般为20ms，道集可以为超道集。其中，该入射角度道集的确定，主要操作为对设定叠加道上目的层段位置，根据目的层段的时窗大小及起、止点的双程反射时间等信息对道集上的各个入射角数据道进行投影，从而得到道集上的相关数据道的目的层段；并对叠加道与各个入射角地震道进行关于目的层的相关系数计算，得到一系列的相关系数数据值；优选相关系数最大的入射角数据道作为叠加道上的相关层位数据的投影道，从而得到相关的某一入射角度的道集数据上的层位数据的投影。原则上，要求最大相关系数数据值大于或等于0.75。其中，相关系数的计算公式如下：
Figure BDA0002979869960000151
式中，Xi及Yi为进行相关系数计算的两种数据的第i个数据值，
Figure BDA0002979869960000161
及
Figure BDA0002979869960000162
分别为两种数据值的等级排序平均值，r的取值范围为0至1。
作为一种优选的技术方案，步骤S3包括以下步骤：
S31,对各个速度分析点上的各个分频道集数据利用相关的层位数据、初始速度谱作为约束，对相关的分频体实施有顺序性的速度谱拾取，得到各个速度分析点上的叠加速度谱；
S32,利用步骤S31获得的各个速度分析点上的叠加速度谱进行叠加处理，得到相关的叠加数据体，再对相关目的层位进行剩余静校正量的计算，并利用剩余静校正量对道集进行校正。
作为一种优选的技术方案，步骤S31中对相关的分频体实施从低频段向高频段的排列顺序的速度谱拾取。
作为一种优选的技术方案，步骤S31中还包括以下步骤：对中低频段、低频段所获得的叠加速度谱、层与层之间相关能量团的聚焦情况等在双程反射时间上进行垂向上速度数据点之间的加密。
作为对步骤S3的详细实施和描述,步骤S3用于以相关速度分析点上的各个分频道集上的层位数据、初始速度谱作为约束，依次对各个速度分析点上各个分频道集数据体进行速度谱的求取，得到各个速度分析点上的叠加速度谱。具体如下：
步骤S3以相关速度分析点上的各个分频道集上的层位数据、初始速度谱进行约束，依次对各个速度分析点上各个分频道集数据体进行精确速度谱的求取，从而得到各个速度分析点上的精确叠加速度谱，包括以下具体步骤：
(3-1)主要根据各个速度分析点上的各个分频道集数据利用相关的层位数据进行约束，依相关的分频体实施有顺序性的速度谱拾取，从而得到各个速度分析点上的相对精确的叠加速度谱。原则上，将第一次拾取的叠加速度谱数据代入第二次叠加速度谱的拾取中，并依相关分频道集数据体从低频段向高频段的排列顺序(本发明所述的顺序性)进行速度谱的拾取并进行更新。一般情况下，在速度谱的拾取中，利用相关的层位数据对其地震反射波进行校平后，所拾取的速度谱点作为相关的速度谱数据点。依次类推，完成各个速度分析点上的精确速度谱点的拾取。其中，速度谱的拾取方法可以采用人工拾取或基于深度强化学习的人工智能拾取方法、将人工与智能拾取方法式进行混合拾取等，具体采用何种方法可以根据实际情况、专家经验及拾取精度等确定。
优选地，对第一次进行精细速度谱拾取的某一分频道集数据进行针对初始速度谱的动校正处理；并根据该速度分析点上常规道集上各个层位数据的偏移距-双程反射时间数据对、坐标信息等对该分频道集数据进行投影，从而得到该分频道集数据的各个层位数据。
优选地，第一次精细速度分析所使用的某一分频道集应该是低频段的分频道集数据，并依低频段、中频段及高频段的频率增高方向进行相关速度分析点上的速度谱拾取。
优选地，鉴于高频段的分频数据体中存在分辨率相对较高的信息，在对高频段的分频道集的速度谱拾取过程中，主要是对中、低频段所获得的叠加速度谱在双程反射时间上进行速度点之间的加密，从而使其相关的反射同相轴达到同相叠加的效果。并根据相关的叠加速度谱的迭代处理，获得更为准确的叠加速度谱。
优选地，如对各个速度分析点上不同层段的分频道集数据，可以对其分别进行叠加速度数据点的拾取，并对不同层段上的同一分频道集上的叠加速度数据点进行合并，从而得到该分频道集所对应的叠加速度谱；并依此类推，完成该速度分析点上的各个分频道集上的叠加速度谱的迭代处理。
(3-2)利用上述步骤获得的各个速度分析点上的叠加速度谱进行叠加处理，得到相关的叠加数据体，再对相关目的层位进行剩余静校正量的计算，并利用剩余静校正量对道集进行校正后，可以再重复步骤1-3，再次进行相关速度分析点上的速度谱的精确拾取，并根据相关叠加速度数据收敛情况终止叠加速度的拾取，从而得到更为精确叠加速度谱。其中，可以按上述相关的处理参数进行处理，也可以根据专家经验、叠加速度数据收敛等情况，确定上述步骤相关的处理参数，进行迭代分析。如可以在迭代中，增加层位的数目、减少分频数据的数目等，从而提高相关的速度分析精度及减少工作量。其中，所述的叠加数据收敛是指经多次叠加速度数据的迭代处理后，相关的叠加数据值趋向某一相对稳定的数据值，并且其相关道集经过叠加处理后的所对应的目的层地震反射振幅的能量值相对较大且稳定。
优选地，如果步骤3-1已经得到最佳速度谱，则步骤3-2可以不实施。对相关最佳速度谱的判断，可以从相关的叠加剖面的反射波连续性、剩余静校正的收敛情况及专家经验等确定。
实施例3
如图1所示，制定工作步骤，对某三维工区的地震资料实施精确速度谱拾取工作，为该研究区地震数据的叠后偏移及叠前深度偏移处理提供相关准确的速度模型。
在步骤S1中确定叠后偏移数据体的各个层位数据，经内插后得到各个速度分析点上的叠加道的层位数据，确定叠前动校正处理后的道集数据体关分频处理参数并经相关带通滤波处理后得到一系列的分频道集数据体。在该步骤中，实际操作中主要为利用测井资料、地质分层资料及叠后偏移资料等，进行井震标定，并且确定相关的目的层。在对本研究区的目的层的分层中，确定14个层位，囊括上、中及下部相关的对地震成像有影响层位数据，并对其进行内插、圆滑等处理后，得到相关的层位数据；并根据相关的速度分析网格，提取网格点上层位数据。此外，对速度分析网格点上提取超级道集数据，设定同一的等间隔的分频处理参数。在实际中，速度分析网格为20线X20道。根据道集的处理频段8hz～70hz，设计相关的分频处理参数，分别为10hz～26hz、24hz～40hz、36hz～52hz、50hz～66hz等四个分频参数，为四个分频数据体。提取速度分析点上的道集数据，并依据四个分频参数应用离散傅里叶变换(DFT)把叠前道集数据变换到频率域，按设计的分频参数利用频率域滤波方法生成一系列分频道集数据体。
在步骤S2中基于道集上的初始速度谱及相关速度分析点上的叠加道上各个层位数据并对相关分频道集数据进行投影，从而得到相关速度分析点上的各个分频道集上的层位数据。在实际操作中，对14个层位数据对各个速度分析点上的道集数据进行层位投影及自动追踪解释，得到相关道集上的14个层位数据；再将14个层位数据按相关信息对四个分频道集数据进行投影，得到四个分频道集上的层位数据。
在步骤S3中以相关速度分析点上的各个分频道集上的层位数据进行约束，依次对各个速度分析点上各个分频道集数据体进行精确速度谱的求取，从而得到各个速度分析点上的精确叠加速度谱。在实际操作中，主要是根据道集上的层位数据的拉平情况，并以此为对速度谱拾取的准确度进行判断；以初始速度谱为背景线，先对10hz～26hz的分频道集数据进行速度谱的拾取，调整相关初始速度谱上的速度数据点在谱中的位置，使其相关层位及层位之间的同相轴进行拉平，从而得到第一次速度谱；再利用第一次速度谱，按相关步骤依次对24hz～40hz、36hz～52hz、50hz～66hz等分频道集进行速度谱的调整及修改。通过对各个速度分析点上的速度谱拾取，从而得到最终的精确速度谱。
利用本发明技术所得到的准确的速度谱成果，与采用常规速度谱分析结果所得到的地层层速度模型(经均方根速度计算所得到)进行相关后续的叠后时间偏移处理后，以两者过井的地震剖面进行对比分析。相关对比成果表明，本发明技术所得到的地震剖面优于常规处理所得到的地震剖面，这也证明了本发明技术对准确速度谱的拾取是有效的。
如上所述，可较好的实现本发明。
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，依据本发明的技术实质，在本发明的精神和原则之内，对以上实施例所作的任何简单的修改、等同替换与改进等，均仍属于本发明技术方案的保护范围之内。

Claims (9)

1.一种拾取精确叠加速度谱的方法，其特征在于，包括以下步骤：
S1,确定相关的叠加道上层位数据,根据分频参数对道集数据进行滤波处理，提取速度分析点上的分频道集数据体；
S2，对速度分析点中不同的分频道集数据体实施各个层位数据的引入，得到相关速度分析点上的各个分频道集上的层位数据；
S3，以相关速度分析点上的各个分频道集上的层位数据、初始速度谱作为约束，依次对各个速度分析点上各个分频道集数据体进行速度谱的迭代性求取，得到各个速度分析点上的叠加速度谱。
2.根据权利要求1或2所述的一种拾取精确叠加速度谱的方法，其特征在于，还包括位于步骤S3之后的以下步骤：
S4，判断各个速度分析网格点的速度谱上叠加速度数据收敛情况，若未收敛，则返回执行步骤S1～S3；若收敛，则结束处理。
3.根据权利要求1或2所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S1包括以下步骤：
S11，准备地质、测井以及地震数据；
S12，建立速度谱拾取约束使用的各个层位数据、初始速度谱；
S13，设定关于叠前道集的分频处理参数。
4.根据权利要求3所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S11中，采用带通滤波器对道集数据进行滤波处理。
5.根据权利要求1或2所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S2包括以下步骤：
S21,基于叠后处理资料中所确定的初始速度谱及速度分析点叠加道上各个层位数据，确定其相关层位的地震反射特征及双程反射时间数据；
S22,对相关速度分析点上的整个频段道集的相关层位数据进行引入，得到该速度分析点道集上的层位数据；
S23，根据速度分析点上的道集上各个层位数据信息提取其分频道集数据，利用道集上的各个层位数据信息对其分频道集数据进行相关层位数据投影，得到相关分频道集的层位数据。
6.根据权利要求5所述的一种拾取精确叠加速度谱的方法，其特征在于，所述的层位数据信息包括但不限于层位的双程反射时间、坐标信息、偏移距信息。
7.根据权利要求1或2所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S3包括以下步骤：
S31,对各个速度分析点上的各个分频道集数据利用相关的层位数据、初始速度谱作为约束，对相关的分频体实施有顺序性的速度谱拾取，得到各个速度分析点上的叠加速度谱；
S32,利用步骤S31获得的各个速度分析点上的叠加速度谱进行叠加处理，得到相关的叠加数据体，再对相关目的层位进行剩余静校正量的计算，并利用剩余静校正量对道集进行校正。
8.根据权利要求7所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S31中对相关的分频体实施从低频段向高频段的排列顺序的速度谱拾取。
9.根据权利要求根据权利要求8所述的一种拾取精确叠加速度谱的方法，其特征在于，步骤S31中还包括以下步骤：对中低频段、低频段所获得的叠加速度谱在双程反射时间上进行垂向上的速度点之间的加密。
CN202110284488.5A 2021-03-17 2021-03-17 一种拾取精确叠加速度谱的方法 Active CN112946742B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110284488.5A CN112946742B (zh) 2021-03-17 2021-03-17 一种拾取精确叠加速度谱的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110284488.5A CN112946742B (zh) 2021-03-17 2021-03-17 一种拾取精确叠加速度谱的方法

Publications (2)

Publication Number Publication Date
CN112946742A true CN112946742A (zh) 2021-06-11
CN112946742B CN112946742B (zh) 2023-10-27

Family

ID=76230296

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110284488.5A Active CN112946742B (zh) 2021-03-17 2021-03-17 一种拾取精确叠加速度谱的方法

Country Status (1)

Country Link
CN (1) CN112946742B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114114399A (zh) * 2021-11-26 2022-03-01 同济大学 一种基于Bayes统计决策的智能速度谱解释与建模方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020042678A1 (en) * 2000-08-07 2002-04-11 Dimitri Bevc Typing picks to horizons in migration velocity analysis
CN104155690A (zh) * 2014-08-08 2014-11-19 中国石油大学(华东) 基于椭球展开的三维地震数据叠加速度求取方法
CN105116446A (zh) * 2015-08-11 2015-12-02 中国石油化工股份有限公司 一种用反射层位约束进行速度分析的方法
CN109143398A (zh) * 2017-06-28 2019-01-04 中国石油化工股份有限公司 一种自动网格层析深度域速度的建模方法
CN109444959A (zh) * 2018-11-01 2019-03-08 科吉思石油技术咨询(北京)有限公司 全频高精度层速度场建立方法
CN109839660A (zh) * 2018-11-08 2019-06-04 成都捷科思石油天然气技术发展有限公司 一种利用叠前道集数据建立速度深度模型的方法
CN111722283A (zh) * 2020-06-28 2020-09-29 成都捷科思石油天然气技术发展有限公司 一种地层速度模型建立方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020042678A1 (en) * 2000-08-07 2002-04-11 Dimitri Bevc Typing picks to horizons in migration velocity analysis
CN104155690A (zh) * 2014-08-08 2014-11-19 中国石油大学(华东) 基于椭球展开的三维地震数据叠加速度求取方法
CN105116446A (zh) * 2015-08-11 2015-12-02 中国石油化工股份有限公司 一种用反射层位约束进行速度分析的方法
CN109143398A (zh) * 2017-06-28 2019-01-04 中国石油化工股份有限公司 一种自动网格层析深度域速度的建模方法
CN109444959A (zh) * 2018-11-01 2019-03-08 科吉思石油技术咨询(北京)有限公司 全频高精度层速度场建立方法
CN109839660A (zh) * 2018-11-08 2019-06-04 成都捷科思石油天然气技术发展有限公司 一种利用叠前道集数据建立速度深度模型的方法
CN111722283A (zh) * 2020-06-28 2020-09-29 成都捷科思石油天然气技术发展有限公司 一种地层速度模型建立方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114114399A (zh) * 2021-11-26 2022-03-01 同济大学 一种基于Bayes统计决策的智能速度谱解释与建模方法
CN114114399B (zh) * 2021-11-26 2023-06-23 同济大学 一种基于Bayes统计决策的智能速度谱解释与建模方法

Also Published As

Publication number Publication date
CN112946742B (zh) 2023-10-27

Similar Documents

Publication Publication Date Title
CN108957549B (zh) 一种辫状河沉积非均质致密砂岩气藏地质建模方法
CN108196305B (zh) 一种山地静校正方法
CN109738945B (zh) 一种利用叠前深度偏移成果直接生成构造图的方法
CN103616723B (zh) 基于avo特征的crp道集真振幅恢复方法
CN111856566B (zh) 湖相滩坝砂体中薄储层预测方法及装置
CN111983683B (zh) 一种少井条件下盐下湖相灰岩储层预测方法和系统
CN106443774A (zh) 一种提高起伏地表叠前深度偏移成像精度的方法
CN109765615A (zh) 一种地层品质因子反演方法及装置
CN106125139B (zh) 一种三维地震数据处理方法及系统
CN110389382A (zh) 一种基于卷积神经网络的油气藏储层表征方法
WO2005026776A1 (en) Wide-offset-range pre-stack depth migration method for seismic exploration
CN105093318B (zh) 一种自适应波动方程波场延拓静校正方法
CN104597485B (zh) 一种微小断层检测方法及断层检测装置
CN103076628A (zh) 一种孔径优化的叠前时间偏移的处理方法
CN112946742A (zh) 一种拾取精确叠加速度谱的方法
CN109884693A (zh) 自适应走向速度谱求取方法和系统
CN110244383B (zh) 基于近地表数据的地质岩性综合模型创建方法
CN111142159B (zh) 中深层异常体地震速度建模方法及系统
CN109143368B (zh) 山前带地震去噪方法
CN112213776B (zh) 叠前道集和vsp资料联合层控q模型建立方法
CN111538080B (zh) 地震成像的方法
CN109839659B (zh) 一种对叠前深度偏移剖面进行迭代优化的方法
CN114428343A (zh) 基于归一化互相关的Marchenko成像方法及系统
CN109581499A (zh) 一种利用各向异性叠前深度偏移生成构造图的方法
CN113917539B (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
CB03 Change of inventor or designer information
CB03 Change of inventor or designer information

Inventor after: Zheng Jian

Inventor after: Qi Xun

Inventor after: Yu Jing

Inventor after: Chen Kelin

Inventor after: He Yanbing

Inventor after: Jing Cui

Inventor after: Nie Zhou

Inventor after: Huang Jun

Inventor after: Luo Hu

Inventor after: Zhang Xiaodan

Inventor before: Zheng Jian

Inventor before: Qi Xun

Inventor before: Yu Jing

Inventor before: He Yanbing

Inventor before: Chen Kelin

Inventor before: Jing Cui

Inventor before: Nie Zhou

Inventor before: Huang Jun

Inventor before: Luo Hu

Inventor before: Zhang Xiaodan

TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20210825

Address after: No.6, building 2-4, No.13, Xinrong Road, Fangcao West 2nd Street, high tech Zone, Chengdu, Sichuan 610000

Applicant after: CHENGDU JIEKESI PETROLEUM NATURAL GAS TECHNOLOGY DEVELOPMENT Co.,Ltd.

Applicant after: SICHUAN CHANGNING NATURAL GAS DEVELOPMENT Co.,Ltd.

Address before: No.6, building 2-4, No.13, Xinrong Road, Fangcao West 2nd Street, high tech Zone, Chengdu, Sichuan 610000

Applicant before: CHENGDU JIEKESI PETROLEUM NATURAL GAS TECHNOLOGY DEVELOPMENT Co.,Ltd.

GR01 Patent grant
GR01 Patent grant