CN113031074A - 一种测井曲线消噪方法、装置、设备及存储介质 - Google Patents

一种测井曲线消噪方法、装置、设备及存储介质 Download PDF

Info

Publication number
CN113031074A
CN113031074A CN202110224658.0A CN202110224658A CN113031074A CN 113031074 A CN113031074 A CN 113031074A CN 202110224658 A CN202110224658 A CN 202110224658A CN 113031074 A CN113031074 A CN 113031074A
Authority
CN
China
Prior art keywords
noise
sequence
logging curve
determining
imf
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
CN202110224658.0A
Other languages
English (en)
Other versions
CN113031074B (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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN202110224658.0A priority Critical patent/CN113031074B/zh
Publication of CN113031074A publication Critical patent/CN113031074A/zh
Application granted granted Critical
Publication of CN113031074B publication Critical patent/CN113031074B/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/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data

Landscapes

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

Abstract

本申请实施例公开了一种测井曲线消噪方法、装置、设备及存储介质，包括：将原始测井曲线数据分解为多个IMF，根据每个IMF确定原始测井曲线的能量密度序列与平均周期序列，绘制原始测井曲线数据对数坐标系下的平均能量‑平均周期散点(下面简称为散点)分布图；根据原始测井曲线构造随机噪音序列，将其分解为多个噪音子序列，将每个噪音子序列分解为多个噪音IMF，确定每个噪音子序列的能量密度序列与平均周期序列，绘制随机噪音序列散点分布图；通过随机噪音序列散点确定噪音分布上下包络线；根据原始测井曲线数据散点分布图与噪音分布上下包络线，确定原始测井曲线数据中噪音成分；根据原始测井曲线数据以及噪音成分，确定消噪后的测井曲线数据。

Description

一种测井曲线消噪方法、装置、设备及存储介质
技术领域
本申请涉及测井技术领域，尤其涉及一种测井曲线消噪方法、装置、设备及存储介质。
背景技术
原始测井曲线(指测井仪器相关传感器采集到的原始地层数据)上会出现许多与地层性质无关的噪音干扰。这些噪音在测井曲线上表现为高频的毛刺干扰。用这些具有毛刺干扰的测井曲线做数字处理，会给计算地质参数带来误差，甚至无法使用。因此，在测井资料预处理中，必须把这些与地层性质无关的毛刺干扰滤除，只保留曲线上反映地层特性的有用信息。
现有的采用滑动平均数字滤波对原始测井曲线进行去噪的方法，需要人为指定滤波参数。而基于HHT的常规去噪方法，针对不同的测井曲线需要进行逐一分析，还需要确定噪音的先验信息。
基于此，需要一种自适应不同原始测井曲线以进行去噪的方法。
发明内容
本申请实施例提供了一种测井曲线消噪方法、装置、设备及存储介质，用于解决如下技术问题：现有对原始测井曲线去噪的方法，难以自适应不同原始测井曲线以进行去噪。
本申请实施例采用下述技术方案：
本申请实施例提供一种测井曲线消噪方法。将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；通过所述随机噪音序列的散点确定噪音分布上下包络线；根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请实施例分别识别原始测井曲线数据的IMF的平均密度和平均周期，以及随机噪音序列的IMF的平均密度和平均周期，判断该阶IMF是否为噪音，进而有效去除噪音。本申请实施例无需指定滤波参数，可自适应识别噪音，适应于不同地区、不同种类的测井曲线。同时，本申请实施例通过散点分布即可识别噪音成分，不依赖使用者的测井曲线数据处理经验，有效提高测井曲线数据的准确性。
在本申请的一种实现方式中，根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据，具体包括：将所述噪音成分对应的IMF阶数进行重构，将重构后的信号作为预测的噪音序列，并确定出所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱；根据所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱，确定所述消噪后的测井曲线数据。
本申请实施例通过构建噪音序列，并依据噪音序列的振幅谱以及原始测井曲线的振幅谱与相位谱，确定消噪后的测井曲线数据。根据测井曲线噪音的分布规律去除噪音，提高测井曲线质量，进而提高测井解释与地质解释质量。
在本申请的一种实现方式中，根据所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱，确定所述消噪后的测井曲线数据，具体包括：将所述原始测井曲线数据的振幅谱减去所述预测的噪音序列的振幅谱，得到消噪的振幅谱；根据所述原始测井曲线数据的相位谱以及所述消噪的振幅谱，构造消噪后的频域信号；对所述消噪后的频域信号进行反傅里叶变换，得到消噪后的测井曲线数据。
在本申请的一种实现方式中，根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分，具体包括：将所述原始测井曲线数据的散点分布图中落在所述噪音分布上下包络线内的散点，确定为所述原始测井曲线数据的噪音成分。
本申请实施例将原始测井曲线数据的散点分布图，与噪音分布上下分布图的位置进行比较，以此确定噪音成分。本申请实施例不仅可以适应不同地区、不同种类的测井曲线，同时也无需指定滤波、无需对先验信息进行确定。因此，本申请实施例不仅简化了去噪的过程，也提高去噪效果。
在本申请的一种实现方式中，确定出IMF的能量密度序列与平均周期序
Figure BDA0002956691290000031
定能量密度序列；其中，En为能量密度，N为当前IMF序列内的点数；Cn(j)为任意一个IMF；根据公式
Figure BDA0002956691290000032
分别确定所述IMF的平均周期；并对所述IMF的平均周期进行对数运算，以确定平均周期序列；其中，T为平均周期；k是IMF的局部极值点总数；Max_index为局部极值点对应的位置序列。
在本申请的一种实现方式中，确定出噪音IMF的能量密度序列与平均周期序列，具体包括：根据公式
Figure BDA0002956691290000033
分别确定噪音IMF的能量密度；根据公式
Figure BDA0002956691290000034
对所述噪音IMF的能量密度进行对数运算，以确定噪音能量密度序列；其中，Noise_En为能量密度，N为当前噪音IMF序列内的点数；Cn(j)为任意一个噪音IMF；根据公式
Figure BDA0002956691290000035
分别确定噪音IMF的平均周期；根据公式
Figure BDA0002956691290000041
对所述噪音IMF的平均周期进行对数运算，确定平均周期序列；其中，Noise_T为所述噪音IMF的平均周期；k是所述噪音IMF的局部极值点总数；Max_index为所述噪音IMF的局部极值点对应的位置序列。
在本申请的一种实现方式中，确定噪音能量密度序列之前，所述方法还包括：根据公式
Figure BDA0002956691290000042
对所述噪音IMF的平均能量值进行校正；其中，
Figure BDA0002956691290000043
为校正后的平均能量值；lnNoise_En为进行对数运算后的平均能量值；
Figure BDA0002956691290000044
为校正系数，其中，所述校正系数是为进行对数运算后的平均能量值的平均值；n为所述噪音IMF的个数；i为所述噪音IMF的序列编号。
本申请实施例提供一种测井曲线消噪装置，包括：确定单元，将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；原始测井曲线散点分布图绘制单元，根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；噪音序列分解单元，根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；噪音序列的散点分布图绘制单元，根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述噪音序列的散点分布图；上下包络线确定单元，通过所述随机噪音序列的散点确定噪音分布上下包络线；噪音成分确定单元，根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；测井曲线数据确定单元，根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请实施例提供一种测井曲线消噪设备，包括：至少一个处理器；以及，与所述至少一个处理器通信连接的存储器；其中，所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够：将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；通过所述随机噪音序列的散点确定噪音分布上下包络线；根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请实施例提供的一种非易失性计算机存储介质，存储有计算机可执行指令，所述计算机可执行指令设置为：将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；通过所述随机噪音序列的散点确定噪音分布上下包络线；根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请实施例采用的上述至少一个技术方案能够达到以下有益效果：本申请实施例分别识别原始测井曲线数据的IMF的平均密度和平均周期，以及随机噪音序列的IMF的平均密度和平均周期，判断该阶IMF是否为噪音，进而有效去除噪音。本申请实施例无需指定滤波参数，可自适应识别噪音，适应于不同地区、不同种类的测井曲线。同时，本申请实施例通过散点分布即可识别噪音成分，不依赖使用者的测井曲线数据处理经验，有效提高测井曲线数据的准确性。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1为本申请实施例提供的一种测井曲线消噪方法流程图；
图2为本申请实施例提供的一种测井曲线消噪整体算法流程图；
图3为本申请实施例提供的一种测井曲线噪音分布规律统计流程图；
图4为本申请实施例提供的一种测井曲线噪音预测流程图；
图5为本申请实施例提供的一种测井曲线噪音消除流程图；
图6为本申请实施例提供的一种渤海某井791m-1012m随钻测井伽马实测数据图；
图7为本申请实施例提供的一种原始测井曲线数据的EMD分解图；
图8为本申请实施例提供的一种噪音周期-能量分布图；
图9为本申请实施例提供的一种预测的噪音图；
图10为本申请实施例提供的一种预测的噪音振幅谱图；
图11为本申请实施例提供的一种原始测井信号振幅谱图；
图12为本申请实施例提供的一种消噪后的振幅谱图；
图13为本申请实施例提供的一种滤波前后对比图；
图14为本申请实施例提供的一种原始信号同滤波信号叠加显示图；
图15为本申请实施例提供的一种测井曲线消噪装置结构示意图；
图16为本申请实施例提供的一种测井曲线消噪设备的结构示意图。
具体实施方式
本申请实施例提供一种测井曲线消噪方法、装置、设备及存储介质。
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本说明书实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。
现有的测井曲线消噪方法，通常采用滑动平均数字滤波对原始测井曲线数据进行去噪。滤波过程中需要人为指定滤波参数，当滤波参数给定过小，滤波效果不明显，毛刺消除效果不好。当滤波参数给定过大，曲线过于平滑，纵向分辨率严重降低。其次，常规滤波方法本质是做低通滤波，但是曲线中的高频信息不仅仅包括毛刺噪音，还可能包含薄地层信息，常规处理方法在滤除毛刺噪音的同时也滤除了有用薄地层信息，降低了曲线的纵向分辨率。
为了解决上述问题，本申请实施例提供了一种测井曲线消噪方法、装置、设备及存储介质。通过原始测井曲线数据的多个IMF的能量密度序列与平均周期序列，绘制散点分布图。以及，通过噪音IMF的能量密度序列与平均周期序列，绘制随机噪音序列的散点分布图与上下包络线。根据散点分布图与上下包络线的位置，确定噪音成分。本申请实施例无需指定滤波参数，可自适应识别噪音，适应于不同地区、不同种类的测井曲线，普及范围较广。
下面通过附图对本申请实施例提出的技术方案进行详细的说明。
图1为本申请实施例提供的一种测井曲线消噪方法流程图。如图1所示，测井曲线消噪方法包括以下步骤：
S101、将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列。
在本申请的一个实施例中，首先通过经验模态分解的方法(Empirical ModeDecomposition，缩写EMD)，对序列长度为N的原始测井曲线数据x(t)进行分解。即将原始测井曲线分解为n个模态分量，得到模态分量集{IMFi,i＝1,2,3,...n}。
具体的，对原始测井曲线利用三次样条插值连接局部极大值，得到上包络线。并且，对原始测井曲线利用三次样条插值连接局部极小值，得到下包络线。求出上、下包络线的均值，得到均值包络线。将原始测井曲线数据减去均值包络线，得到中间信号。在中间信号满足分解停止规则的情况下，该中间序列就是一个IMF。如果不满足，重复上述过程。
在本申请的一个实施例中，分解停止规则定义为两次操作的归一化平方差：
Figure BDA0002956691290000081
其中，SDk为筛分门限值，在SDk小于预设值的情况下，分解过程停止。hk(t)为当前正在计算的IMF分量，hk-1为前一次计算的IMF分量。
在本申请的一个实施例中，将得到第一个IMF，记为IMF1,用原始测井曲线数据减去IMF1，将得到的结果作为新的原始信号。再通过确定上、下包络线，以及确定均值包络线与中间信号的方法，得到第二个IMF，以此类推，完成EMD分解。当剩余项变得足够小或变成单调函数，从它之中不能再分解任何IMF，即可完成原始测井曲线数据的分解。最终得到：
Figure BDA0002956691290000082
需要说明的是，x(t)为原始测井曲线数据，n为经验模态的数量，IMFj为模态分量，rn为剩余项。
例如，图6为本申请实施例提供的一种渤海某井791m-1012m随钻测井伽马实测数据图。如图所示，横坐标为测井深度，纵坐标为测井曲线数据的幅值。由图6可以明显看出，测井曲线数据中带有部分噪音信号。
图7为本申请实施例提供的一种原始测井曲线数据的EMD分解图。如图所示，将原始测井曲线分解为8个模态分量IMF1-IMF8，最后一项为剩余项。其中每一个模态分量图中的横坐标为测井深度，纵坐标分别为每一个模态分量信号的幅值。
需要说明的是，通过EMD分解，能使复杂信号分解为有限个本征模函数(IntrinsicMode Function，简称IMF)，所分解出来的各IMF分量包含了原信号的不同时间尺度的局部特征信号,特别适用于非线性非平稳信号的分析处理。
在本申请的一个实施例中，对于每一个模态分量IMFi，求解能量密度Ei以及平均周期Ti，构成点集{(Ei,Ti)，i＝1,2,3...,n}。
具体的，根据公式
Figure BDA0002956691290000091
分别确定每一个IMF的能量密度。根据公式
Figure BDA0002956691290000092
对IMF的能量密度进行对数运算，以确定能量密度序列。其中，En为能量密度，N为当前IMF序列内的点数，Cn(j)为任意一个IMF。
具体的，根据公式
Figure BDA0002956691290000093
分别确定IMF的平均周期，并对IMF的平均周期进行对数运算，以确定平均周期序列。其中，T为平均周期，k是IMF的局部极值点总数，Max_index为局部极值点对应的位置序列。
在本申请的一个实施例中，通过对IMF依次计算能量密度对数，得到能量密度序列{ln Ei,i＝0,1,2,3...n}。通过对IMF依次计算平均周期对数，得到平均周期序列{ln Ti,i＝0,1,2,3...n}。
S102、根据IMF的能量密度序列与平均周期序列，绘制原始测井曲线数据的散点分布图。
在本申请的一个实施例中，根据原始测井曲线数据分解后，得到的多个IMF的能量密度序列与平均周期序列，绘出平均周期和能量密度构成的散点{(lnTi，lnEi),i＝0,1,2,3...n}的空间分布图。
S103、根据原始测井曲线的长度构造长随机噪音序列，将长随机噪音序列分解为多个噪音子序列，每个子序列长度和原始测井曲线序列长度相等，将各噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列。
在本申请的一个实施例中，根据原始测井曲线信号长度，生成长随机噪音序列n(t)，序列长度为K×N，例如，可以使用python语言中的随机函数随机生成长噪音序列。其中，K为子序列的数量，N为信号长度。
在本申请的一个实施例中，对长随机噪音序列n(t)进行二范数正则化处理。将处理后的长随机噪音序列进行序列分解，得到K个子序列。其中，每个子序列长度与原始测井曲线信号长度相同。
在本申请的一个实施例中，将每个噪音子序列进行EMD分解，确定出多个噪音IMF。所有噪音子序列分解得到的IMF可以记为{Noise_imfs_i,i＝0,1,2,3....}。
在本申请的一个实施例中，对于Nosie_imf_i序列中的每一个IMF，计算其平均周期和能量密度，并确定出噪音IMF的能量密度序列与平均周期序列。
具体的，根据公式
Figure BDA0002956691290000101
分别确定噪音IMF的能量密度。
根据公式
Figure BDA0002956691290000102
对噪音IMF的能量密度进行对数运算，以确定噪音能量密度序列。
其中，Noise_En为能量密度，N为当前噪音IMF序列内的点数，Cn(j)为任意一个噪音IMF。
在本申请的一个实施例中，由于地层随机统计均值不为零。因此需要对上述噪音的平均能量值进行校正，以平均能量值平均数为校正系数，将上述噪音所有点的平均能量值加上校正系数。
具体的，根据公式
Figure BDA0002956691290000111
对所述噪音IMF的平均能量值进行校正。
其中，公式
Figure BDA0002956691290000112
为校正后的平均能量值。
公式
Figure BDA0002956691290000113
为进行对数运算后的平均能量值。
公式
lnNoise_En
为校正系数。其中，校正系数是为进行对数运算后的平均能量值的平均值；n为噪音IMF的个数；i为噪音IMF的位置序列编号。
在本申请的一个实施例中，根据公式
Figure BDA0002956691290000114
分别确定噪音IMF的平均周期。得到噪音IMF的平均周期后，根据公式
Figure BDA0002956691290000115
对噪音IMF的平均周期进行对数运算，确定平均周期序列。其中，Noise_T为噪音IMF的平均周期，k是噪音IMF的局部极值点总数。Max_index为噪音IMF的局部极值点对应的位置序列。
在本申请的一个实施例中，对噪音IMF依次进行能量密度值计算，并对噪音能量密度值进行对数运算，可以得到噪音能量密度序列
Figure BDA0002956691290000116
对噪音IMF依次计算平均周期，并对噪音平均周期进行对数运算，得到噪音平均周期序列{lnNoise_Ti,i＝0,1,2,3...n}。
S104、根据噪音IMF的能量密度序列与平均周期序列，绘制随机噪音序列的散点分布图。
在本申请的一个实施例中，统计得到的随机噪音序列的能量密度序列与平均周期序列，将序列点
Figure BDA0002956691290000121
绘制成图。
图8为本申请实施例提供的一种噪音周期-能量分布图，如图8所示，点A、点B、点C、点D、点E、点F、点G、点H分别为原始测井曲线数据，经处理后生成的信号分量。表明原始信号在lnT-lnE坐标系下的分布规律。图8中除了原始测井信号生成的信号分量之外，其它的点为构造的超长的噪音序列，经处理后生成的噪音分量。表明随机噪音在lnT-lnE坐标系下的分布规律。
在本申请的一个实施例中，如果原始测井曲线数据，经处理后生成的信号分量落在了噪音分量点集中的区域，代表该分量是噪音，如果信号分量点偏离噪音分量点区域，代表该分量是有用信号。
S105、通过随机噪音序列的散点确定噪音分布上下包络线。
如图8所示，图中上下两条包络线分别为噪音分布的上下包络线。本申请实施例可以利用三次样条插值连接局部极大值，得到噪音分量的上包络线。以及利用三次样条插值连接局部极小值，得到噪音分量的下包络线条。
S106、根据原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出原始测井曲线数据中的噪音成分。
在本申请的一个实施例中，将原始测井曲线数据的散点分布图中落在噪音分布上下包络线内的散点，确定为原始测井曲线数据的噪音成分。将落在噪音分布边界外的信号点，认为是有用成分。
需要说明的是，原始测井曲线数据的散点分布图中落在上下包络线上的散点，属于临界状态，可能包含噪音。因此可以根据测量需要对临界状态的散点进行判定。
S107、根据原始测井曲线数据以及噪音成分，确定出消噪后的测井曲线数据。
在本申请的一个实施例中，将噪音成分对应的IMF阶数进行重构，将重构后的信号作为预测的噪音序列y(t)。
图9为本申请实施例提供的一种预测的噪音序列图，如图9所示，横坐标为测井深度，纵坐标为噪音序列的振幅。通过图9可以明显看出噪音序列的振幅变化趋势。
在本申请的一个实施例中，对预测的噪音序列进行傅里叶变换，确定出预测的噪音序列的振幅谱|Y(ω)|。以及对原始信号进行傅里叶变换，得到原始测井曲线数据的振幅谱|X(ω)|与相位谱
Figure BDA0002956691290000131
图10为本申请实施例提供的一种预测的噪音振幅谱图。如图10所示，横坐标为预测的噪音的频率，纵坐标为预测的噪音的振幅。通过图10可以明显看出预测的噪音振幅谱的变化趋势。图11为本申请实施例提供的一种原始测井信号振幅谱图。如图11所示，横坐标为原始测井信号的频率，纵坐标为原始测井信号的振幅。通过图11可以明显看出原始测井信号振幅谱的变化趋势。
在本申请的一个实施例中，根据预测的噪音序列的振幅谱，以及原始测井曲线数据的振幅谱与相位谱，确定消噪后的测井曲线数据。
具体的，将原始测井曲线数据的振幅谱减去预测的噪音序列的振幅谱，得到消噪的振幅谱|S(ω)|。根据原始测井曲线数据的相位谱
Figure BDA0002956691290000132
以及消噪的振幅谱|S(ω)|，构造消噪后的频域信号。根据函数
Figure BDA0002956691290000133
对消噪后的频域信号进行反傅里叶变换，得到消噪后的测井曲线数据
Figure BDA0002956691290000134
图12为本申请实施例提供的一种消噪后的振幅谱图。如图12所示，横坐标为消噪后信号的频率，纵坐标为消噪后信号的振幅。图13为本申请实施例提供的一种滤波前后对比图。如图13所示，上图为原始测井曲线信号图，其横坐标为测井深度，纵坐标为振幅。下图为去噪后的信号图，其横坐标为测井深度，纵坐标为振幅。图14为本申请实施例提供的一种原始信号同滤波信号叠加显示图。如图14所示，横坐标为测井深度，纵坐标为信号振幅。图中颜色较浅的波形为原始信号波形，叠加在原始信号波形上面的颜色较深的波形，为去噪后的测井曲线信号。
本申请实施例分别识别原始测井曲线数据的IMF的平均密度和平均周期，以及噪音序列的IMF的平均密度和平均周期，判断该阶IMF是否为噪音，进而有效去除噪音。本申请实施例无需指定滤波参数，可自适应识别噪音，适应于不同地区、不同种类的测井曲线。通过散点分布即可识别噪音成分，不依赖使用者的测井曲线数据处理经验，有效提高测井曲线数据的准确性。
图2为本申请实施例提供的一种测井曲线消噪整体算法流程图，如图2所示，测井曲线消噪整体算法流程如下：
在本申请的一个实施例中，对原始测井曲线数据进行EMD分解，得到IMF序列。根据IMF的能量密度序列与平均周期序列，绘制原始测井曲线数据的散点分布图。
在本申请的一个实施例中，噪音预测模块根据原始测井曲线数据构造噪音序列。将随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF。根据噪音IMF的能量密度序列与平均周期序列，绘制噪音序列的散点分布图。
在本申请的一个实施例中，噪音消除模块根据原始测井曲线数据以及噪音成分，对原始测井曲线数据进行消噪，得到消噪后的测井曲线数据。根据消噪后的测井曲线数据，重构测井曲线。
图3为本申请实施例提供的一种测井曲线噪音分布规律统计流程图，如图3所示，测井曲线噪音分布规律统计流程如下：
在本申请的一个实施例中，分别确定噪音序列的平均周期与平均能量，并对平均能量进行能量校正。
在本申请的一个实施例中，根据确定出的噪音序列的平均周期与校正后的能量，分别对平均周期与校正后的能量进行对数计算，根据对数计算的结果进行对数坐标系绘图。
在本申请的一个实施例中，在对数坐标系中绘制散点分布上包络线与下包络线。
在本申请的一个实施例中，根据绘制出的对数坐标系以及上下包络线，可以得到噪音各模态分量的平均周期-能量空间分布规律。
图4为本申请实施例提供的一种测井曲线噪音预测流程图，如图4所示，测井曲线噪音预测流程如下：
在本申请的一个实施例中，将原始测井曲线进行EMD分解，得到IMF序列。确定出IMF的能量密度序列与平均周期序列，根据IMF的能量密度序列与平均周期序列，绘制原始测井曲线数据的散点分布图。
在本申请的一个实施例中，根据原始测井曲线构造噪音序列，将随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF序列，并确定出噪音IMF的能量密度序列与平均周期序列。根据噪音IMF的能量密度序列与平均周期序列，绘制噪音序列的散点分布图。
在本申请的一个实施例中，判断原始测井曲线某阶模态分量是否落在噪音包络线内。在原始测井曲线数据的散点分布图中落在噪音分布上下包络线内的散点，确定为原始测井曲线数据的噪音成分。在原始测井曲线数据的散点分布图中落在噪音分布上下包络线外的散点，确定为原始测井曲线数据的有用信号。
在本申请的一个实施例中，记录下为噪音的模态分量阶数，将噪音成分对应模态分量阶数进行重构，将重构后的信号作为预测的噪音序列。
图5为本申请实施例提供的一种测井曲线噪音消除流程图，如图5所示。测井曲线噪音消除流程如下：
在本申请的一个实施例中，将噪音成分对应的IMF阶数进行重构，将重构后的信号作为预测的噪音序列，并确定出预测的噪音序列的振幅谱。同时根据原始信号，确定原始信号的振幅谱与相位谱。
在本申请的一个实施例中，将原始测井曲线数据的振幅谱减去预测的噪音序列的振幅谱，得到消噪的振幅谱。
在本申请的一个实施例中，利用原始信号的相位谱
Figure BDA0002956691290000151
和消噪后的振幅谱|S(ω)|构造消噪后的频域信号，对该频域信号进行反傅里叶变换，得到消噪后的时域信号。
图15为本申请实施例提供的一种测井曲线消噪装置结构示意图，装置包括：确定单元S301、原始测井曲线散点分布图绘制单元S302、噪音序列分解单元S303、噪音序列的散点分布图绘制单元S304、上下包络线确定单元S305、噪音成分确定单元S306、测井曲线数据确定单元S307。
确定单元S301，将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列。
原始测井曲线散点分布图绘制单元S302，根据IMF的能量密度序列与平均周期序列，绘制原始测井曲线数据的散点分布图。
噪音序列分解单元S303，根据原始测井曲线构造随机噪音序列，将随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列。
噪音序列的散点分布图绘制单元S304，根据噪音IMF的能量密度序列与平均周期序列，绘制噪音序列的散点分布图。
上下包络线确定单元S305，通过随机噪音序列的散点确定噪音分布上下包络线。
噪音成分确定单元S306，根据原始测井曲线数据的散点分布图与噪音分布上下包络线，确定出原始测井曲线数据中的噪音成分。
测井曲线数据确定单元S307，根据原始测井曲线数据以及噪音成分，确定出消噪后的测井曲线数据。
图16为本申请实施例提供的一种测井曲线消噪设备的结构示意图。
本申请实施例提供的一种测井曲线消噪设备，包括：
至少一个处理器；以及，
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够：
将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列。
根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图。
根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列。
根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图。
通过所述随机噪音序列的散点确定噪音分布上下包络线。
根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分。
根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请实施例提供一种非易失性计算机存储介质，存储有计算机可执行指令，所述计算机可执行指令设置为：
将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列。
根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图。
根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列。
根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图。
通过所述随机噪音序列的散点确定噪音分布上下包络线。
根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分。
根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
本申请中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置、设备、非易失性计算机存储介质实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
上述对本申请特定实施例进行了描述。其它实施例在所附权利要求书的范围内。在一些情况下，在权利要求书中记载的动作或步骤可以按照不同于实施例中的顺序来执行并且仍然可以实现期望的结果。另外，在附图中描绘的过程不一定要求示出的特定顺序或者连续顺序才能实现期望的结果。在某些实施方式中，多任务处理和并行处理也是可以的或者可能是有利的。
以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请的实施例可以有各种更改和变化。凡在本申请实施例的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。

Claims (10)

1.一种测井曲线消噪方法，其特征在于，所述方法包括：
将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；
根据所述能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；
根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；
根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；
通过所述随机噪音序列的散点确定噪音分布上下包络线；
根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；
根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
2.根据权利要求1所述的一种测井曲线消噪方法，其特征在于，所述根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据，具体包括：
将所述噪音成分对应的IMF阶数进行重构，将重构后的信号作为预测的噪音序列；
确定出所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱；
根据所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱，确定所述消噪后的测井曲线数据。
3.根据权利要求2所述的一种测井曲线消噪方法，其特征在于，所述根据所述预测的噪音序列的振幅谱，以及所述原始测井曲线数据的振幅谱与相位谱，确定所述消噪后的测井曲线数据，具体包括：
将所述原始测井曲线数据的振幅谱减去所述预测的噪音序列的振幅谱，得到消噪的振幅谱；
根据所述原始测井曲线数据的相位谱以及所述消噪的振幅谱，构造消噪后的频域信号；
对所述消噪后的频域信号进行反傅里叶变换，得到消噪后的测井曲线数据。
4.根据权利要求1所述的一种测井曲线消噪方法，其特征在于，所述根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分，具体包括：
将所述原始测井曲线数据的散点分布图中落在所述噪音分布上下包络线内的散点，确定为所述原始测井曲线数据的噪音成分。
5.根据权利要求1所述的一种测井曲线消噪方法，其特征在于，所述确定出IMF的能量密度序列与平均周期序列，具体包括：
根据公式
Figure FDA0002956691280000021
分别确定所述IMF的能量密度；
根据公式
Figure FDA0002956691280000022
对所述IMF的能量密度进行对数运算，以确定能量密度序列；其中，En为能量密度，N为当前IMF序列内的点数；Cn(j)为任意一个IMF；
根据公式
Figure FDA0002956691280000023
分别确定所述IMF的平均周期；并对所述IMF的平均周期进行对数运算，以确定平均周期序列；其中，T为平均周期；k是IMF的局部极值点总数；Max_index为局部极值点对应的位置序列。
6.根据权利要求1所述的一种测井曲线消噪方法，其特征在于，所述确定出噪音IMF的能量密度序列与平均周期序列，具体包括：
根据公式
Figure FDA0002956691280000031
分别确定噪音IMF的能量密度；
根据公式
Figure FDA0002956691280000032
对所述噪音IMF的能量密度进行对数运算，以确定噪音能量密度序列；其中，Noise_En为能量密度，N为当前噪音IMF序列内的点数；Cn(j)为任意一个噪音IMF；
根据公式
Figure FDA0002956691280000033
分别确定噪音IMF的平均周期；
根据公式
Figure FDA0002956691280000034
对所述噪音IMF的平均周期进行对数运算，确定平均周期序列；其中，Noise_T为所述噪音IMF的平均周期；k是所述噪音IMF的局部极值点总数；Max_index为所述噪音IMF的局部极值点对应的位置序列。
7.根据权利要求6所述的一种测井曲线消噪方法，其特征在于，所述确定噪音能量密度序列之前，所述方法还包括：
根据公式
Figure FDA0002956691280000035
对所述噪音IMF的平均能量值进行校正；其中，
Figure FDA0002956691280000036
为校正后的平均能量值；lnNoise_En为进行对数运算后的平均能量值；
Figure FDA0002956691280000037
为校正系数，其中，所述校正系数为进行对数运算后的平均能量值的平均值；n为所述噪音IMF的个数；i为所述噪音IMF的序列编号。
8.一种测井曲线消噪装置，包括：
确定单元，将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；
原始测井曲线散点分布图绘制单元，根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；
噪音序列分解单元，根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；
噪音序列的散点分布图绘制单元，根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述噪音序列的散点分布图；
上下包络线确定单元，通过所述随机噪音序列的散点确定噪音分布上下包络线；
噪音成分确定单元，根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；
测井曲线数据确定单元，根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
9.一种测井曲线消噪设备，包括：
至少一个处理器；以及，
与所述至少一个处理器通信连接的存储器；其中，
所述存储器存储有可被所述至少一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器能够：
将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；
根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；
根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；
根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；
通过所述随机噪音序列的散点确定噪音分布上下包络线；
根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；
根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
10.一种非易失性计算机存储介质，存储有计算机可执行指令，所述计算机可执行指令设置为：
将原始测井曲线数据分解为多个IMF，并确定出IMF的能量密度序列与平均周期序列；
根据所述IMF的能量密度序列与平均周期序列，绘制所述原始测井曲线数据的散点分布图；
根据所述原始测井曲线构造随机噪音序列，将所述随机噪音序列分解为多个噪音子序列，将每个噪音子序列进行EMD分解，确定出多个噪音IMF，并确定出噪音IMF的能量密度序列与平均周期序列；
根据所述噪音IMF的能量密度序列与平均周期序列，绘制所述随机噪音序列的散点分布图；
通过所述随机噪音序列的散点确定噪音分布上下包络线；
根据所述原始测井曲线数据的散点分布图与所述噪音分布上下包络线，确定出所述原始测井曲线数据中的噪音成分；
根据所述原始测井曲线数据以及所述噪音成分，确定出消噪后的测井曲线数据。
CN202110224658.0A 2021-03-01 2021-03-01 一种测井曲线消噪方法、装置、设备及存储介质 Active CN113031074B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110224658.0A CN113031074B (zh) 2021-03-01 2021-03-01 一种测井曲线消噪方法、装置、设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110224658.0A CN113031074B (zh) 2021-03-01 2021-03-01 一种测井曲线消噪方法、装置、设备及存储介质

Publications (2)

Publication Number Publication Date
CN113031074A true CN113031074A (zh) 2021-06-25
CN113031074B CN113031074B (zh) 2021-12-07

Family

ID=76466328

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110224658.0A Active CN113031074B (zh) 2021-03-01 2021-03-01 一种测井曲线消噪方法、装置、设备及存储介质

Country Status (1)

Country Link
CN (1) CN113031074B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102998706A (zh) * 2012-11-23 2013-03-27 中国石油大学(北京) 一种衰减地震数据随机噪声的方法及系统
CN108999606A (zh) * 2018-08-10 2018-12-14 常州工学院 一种随钻方位电磁波电阻率测井曲线降噪方法
CN109061736A (zh) * 2018-08-09 2018-12-21 中国石油天然气股份有限公司 一种测井资料周期性噪音消除方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102998706A (zh) * 2012-11-23 2013-03-27 中国石油大学(北京) 一种衰减地震数据随机噪声的方法及系统
CN109061736A (zh) * 2018-08-09 2018-12-21 中国石油天然气股份有限公司 一种测井资料周期性噪音消除方法及系统
CN108999606A (zh) * 2018-08-10 2018-12-14 常州工学院 一种随钻方位电磁波电阻率测井曲线降噪方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
ZHAOHUA WU 等: "A study of the characteristics of white noise using the empirical mode decomposition method", 《THE ROYAL SOCIETY》 *
刘建建: "Hilbert-Huang变换在测井资料处理中的应用研究", 《中国优秀博硕士学位论文全文数据库（硕士）·基础科学辑》 *
高云超 等: "基于经验模式分解的自适应去噪算法", 《计算机工程与应用》 *

Also Published As

Publication number Publication date
CN113031074B (zh) 2021-12-07

Similar Documents

Publication Publication Date Title
CN102998706B (zh) 一种衰减地震数据随机噪声的方法及系统
CN108416740B (zh) 一种消除椒盐噪声的迭代自适应中值滤波方法
CN110648290A (zh) 一种基于sure参数优化的双核非局部均值图像去噪方法
CN112446323B (zh) 基于改进emd模态混叠和端点效应的hht谐波分析方法
CN102930149B (zh) 基于pca和emd的传感器网络感知信息去噪处理方法
CN111007566B (zh) 一种曲率驱动扩散全卷积网络地震数据坏道重建与去噪方法
CN112508810A (zh) 非局部均值盲图像去噪方法、系统及装置
CN110490814B (zh) 基于光滑秩约束的混合噪声去除方法、系统及存储介质
CN108470018A (zh) 基于经验模态分解分解的内蕴模式函数的光滑方法及装置
CN106680876A (zh) 一种地震数据联合去噪方法
CN113723171B (zh) 基于残差生成对抗网络的脑电信号去噪方法
CN115700544A (zh) 一种联合经验模态分解及小波软阈值的色谱信号去噪方法
CN117454085B (zh) 一种车辆在线监控方法及系统
CN113031074B (zh) 一种测井曲线消噪方法、装置、设备及存储介质
CN114077852A (zh) 一种强噪声光谱信号的智能去噪方法
CN110703089B (zh) 一种用于低频振荡Prony分析的小波阈值去噪方法
CN113375065B (zh) 管道泄漏监测中趋势信号的消除方法及装置
CN107907542B (zh) 一种ivmd及能量估计相结合的dspi相位滤波方法
CN115859054A (zh) 基于mic和ceemdan的水电机组尾水管压力脉动数据滤波方法
CN116148935A (zh) 一种基于自适应自编码器的磁共振随机噪声抑制方法
CN115859048A (zh) 一种局放信号的噪声处理方法及装置
CN115795272A (zh) 一种基于分数阶迭代离散小波变换的谱线去噪方法
CN115497492A (zh) 一种基于全卷积神经网络的实时语音增强方法
CN114881883A (zh) 一种红外图像多维降噪方法、存储介质及装置
CN113742985B (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
GR01 Patent grant
GR01 Patent grant