CN115170428A - 一种声波远探测成像图的降噪方法 - Google Patents
一种声波远探测成像图的降噪方法 Download PDFInfo
- Publication number
- CN115170428A CN115170428A CN202210845762.6A CN202210845762A CN115170428A CN 115170428 A CN115170428 A CN 115170428A CN 202210845762 A CN202210845762 A CN 202210845762A CN 115170428 A CN115170428 A CN 115170428A
- Authority
- CN
- China
- Prior art keywords
- noise
- imaging
- wave
- data
- field
- 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.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 92
- 238000000034 method Methods 0.000 title claims abstract description 38
- 230000009467 reduction Effects 0.000 title claims abstract description 31
- 238000001514 detection method Methods 0.000 title claims abstract description 15
- 238000012549 training Methods 0.000 claims abstract description 23
- 238000012545 processing Methods 0.000 claims abstract description 21
- 238000013528 artificial neural network Methods 0.000 claims abstract description 9
- 238000007781 pre-processing Methods 0.000 claims abstract description 7
- 238000003062 neural network model Methods 0.000 claims abstract description 4
- 230000001902 propagating effect Effects 0.000 claims description 12
- 238000005070 sampling Methods 0.000 claims description 12
- 238000004613 tight binding model Methods 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 238000010606 normalization Methods 0.000 claims description 10
- 230000002441 reversible effect Effects 0.000 claims description 10
- 238000012795 verification Methods 0.000 claims description 10
- 230000004913 activation Effects 0.000 claims description 7
- 238000004088 simulation Methods 0.000 claims description 7
- 238000012360 testing method Methods 0.000 claims description 7
- 230000005284 excitation Effects 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 238000013527 convolutional neural network Methods 0.000 claims description 5
- 238000005520 cutting process Methods 0.000 claims description 4
- 238000011478 gradient descent method Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000012935 Averaging Methods 0.000 claims description 3
- 150000001875 compounds Chemical class 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000004044 response Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000011065 in-situ storage Methods 0.000 claims 1
- 230000000694 effects Effects 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 3
- 230000005012 migration Effects 0.000 description 3
- 238000013508 migration Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000010200 validation analysis Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000010219 correlation analysis Methods 0.000 description 1
- 238000013434 data augmentation Methods 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 238000010859 live-cell imaging Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20081—Training; Learning
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20084—Artificial neural networks [ANN]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30168—Image quality inspection
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Molecular Biology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- Artificial Intelligence (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种声波远探测成像图的降噪方法，所述该方法包括以下步骤：步骤1：准备数据集；步骤2：构建降噪网络模型；步骤3：数据预处理；步骤4：训练网络模型；步骤5：现场成像图处理。本发明通过图像处理的角度对声波远探测成像图进行直接降噪，构建降噪神经网络模型，训练端到端的降噪神经网络，不需要复杂的参数设置，在保留成像图边缘信息的同时，自动快速地移除成像图中干扰噪声，提高了成像图的质量。
Description
技术领域
本发明属于应用地球物理声学测井技术领域，具体涉及一种声波远探测成像图的降噪方法。
背景技术
声波远探测测井技术在井中利用声源向井外地层中辐射弹性波，经井外非均匀体反射后被井中接收器采集，之后通过偏移成像算法得到井外地质结构成像图。
在弹性波传播的过程中，相比于沿井壁传播的纵波、横波、斯通利波等直达模式波，有用的反射波由于长距离的几何和黏弹衰减，其幅度极低，甚至低于井下复杂环境引发的背景噪声。反射信号的低信噪比问题严重降低了目前的声波远探测成像效果，在成像图中引发的强烈噪声给有效储层的识别引入了巨大的不确定性，而错误的储层评估将会导致钻井轨迹严重脱离靶区，造成巨大的经济损失。
为了克服该问题，基于双数复小波变换和慢度时间相关法的多尺度相关分析方法综合利用尺度-时间-慢度域信息分离反射波和滑行波(Tao et al.,2008；Wang et al.,2011)，但该方法受到所选择的小波和参数的限制。参数估计法也被用于反射波提取(Tang,2007；Wang et al.,2011；Cai et al.,2020)，其基于设定的物理模型，利用最小二乘方法对直达波场频谱进行估计，然后从总波场中进行扣除，但该方法的效果严重依赖于物理模型的准确性。武宏亮等将偶极横波远探测中的干扰噪声分为相干噪声和非相干噪声(Wu etal.,2021)，在对噪声产生机理和典型特征研究的基础上分别采用对应的处理方法进行压制，提高了反射横波的提取精度。然而在实际测量情况下反射波幅度微弱，且常与井中直达模式波混叠，使得反射波场分离变得异常困难，导致经过反射波提取算法处理后的波场中仍然存有大量残余信号，其幅度甚至高于提取的反射波。这部分残余信号经过后续的偏移成像算法处理后，将会在成像图中呈现为非规则的背景噪声，降低成像图质量，误导后续的储层识别。
发明内容
本发明针对上述的不足之处提供一种能够智能准确地移除成像图中干扰噪声的声波远探测成像图的降噪方法。
本发明目的是这样实现的：一种声波远探测成像图的降噪方法，其特征在于：所述该方法包括以下步骤：
步骤1：准备数据集；
步骤2：构建降噪网络模型；
步骤3：数据预处理；
步骤4：训练网络模型；
步骤5：现场成像图处理。
优选的，所述步骤1中准备数据集的具体操作为：
步骤1-1：利用均匀地层作为背景介质，设定其尺寸、纵横波速度及密度值，设定非均匀介质的纵横波速度和密度值，并利用随机函数确定其尺寸、形状、位置，将非均匀区域叠加在均匀地层之上，构建地质模型；
步骤1-2：通过求解均匀介质的波动方程评价入射波场，然后对非均匀区域进行离散化，对各网格声场响应进行叠加评价散射波场，通过时域积分方程法评价所有时刻的接收波形数据；
步骤1-3：利用偏移成像算法将中得到的时域波形向倾斜同相轴下倾方向归位，得到井周地质结构的标签成像图；
步骤1-4：给标签成像图随机添加不同等级的白高斯噪声，得到噪声成像图；
步骤1-5：重复上步骤1-1、1-2、1-3、1-4，得到大量成对的标签和噪声成像图，建立大规模的理论模拟数据集。
优选的，所述步骤2中构建降噪网络模型采用深层卷积神经网络，构建深层卷积神经网络包括以下步骤：
步骤2-1：通过间隔取值的方式的对尺寸为M×N×1的输入噪声成像图进行可逆降采样，生成尺寸为
的含噪声子图，其中最后一维表示通道；
步骤2-2：建立与含噪声子图相同尺寸
的二维噪声等级图，并全部填充为所添加的白高斯噪声的等级值，将含噪声子图与噪声等级图进行拼接，得到尺寸为
的网络输入；
步骤2-3：利用尺寸为3×3、步长为1×1、维度为5的卷积核对网络输入进行处理，并采用ReLU激活函数进行非线性变换，然后通过尺寸为3×3、步长为1×1、维度为64的卷积核、ReLU激活函数和批归一化层构建卷积层，并连续累加14组卷积层在同一分辨率进行特征提取，充分保存图像中的位置信息；利用尺寸为3×3、步长为1×1、维度为4的卷积核将提取特征转化为降噪后的子图，通过可逆降采样的反操作合并各子图得到降噪后成像图。
优选的，所述步骤3中数据预处理的具体操作为：
步骤3-1：按照75％、15％、10％的比例将数据集分解为训练集、验证集和测试集；
步骤3-2：对图像进行归一化：
其中μ表示数据的均值，σ为数据的标准差；
步骤3-3：随机选择偏移、裁剪、缩放、翻转的方式对数据进行增广。
优选的，所述步骤4中训练网络模型的具体操作为：
步骤4-1：按照设定的批次采样数从训练集中进行随机采样，将噪声成像图和噪声等级值输入给步骤4中构建的降噪神经网络模型；
步骤4-2：联合网络输出与对应的标签数据计算L1损失函数：
其中N为批次采样值，p为网络输出，y为标签成像图，|·|为绝对值运算符；
步骤4-3：利用随机梯度下降法将计算的L1损失函数反向传播，更新神经网络参数；
步骤4-4：重复步骤4-1至步骤4-3对网络进行连续优化，直至训练集中所有数据均被采样完成，对验证集数据进行处理，计算L1损失函数与信噪比指数：
SNR＝10log10(‖p‖2/‖p-y‖2)；
其中，SNR为信噪比指数；
步骤4-5：重复步骤4-4，直至迭代次数达到200，或者5次迭代过程中验证集上的信噪比指数上升小于0.01；
步骤4-6：将训练好的网络应用于测试数据集，评价网络的泛化性能。
优选的，所述步骤5中现场成像图处理的具体操作为：
步骤5-1：将步骤4中训练结束的模型参数进行固化和保存；
步骤5-2：首先按照理论模拟成像图的尺寸对成像图进行滑动切片，在边界处存在重叠区域；
步骤5-3：将分别归一化后的切片与设定的噪声等级值输入给降噪网络，将各输出进行反归一化后重新拼接，对于重叠区域则采用取平均值的方式，最终得到降噪后的成像图。
优选的，所述步骤1-2中通过求解均匀介质的波动方程评价入射波场的具体操作为：
步骤1-2-1：当井外地层中存在传播声速变化的非均匀介质时，井中声源激发的弹性波在该区域内将会改变其原始传播方向，产生散射声场；
步骤1-2-2：在弱散射条件下，非均匀介质速度变化引起的局部不均匀性，散射波的传播方程由入射波场与介质的非均匀性确定：
其中u为标量位移势，V0为非均匀介质的背景速度；
步骤1-2-3：基于Born近似原理，总声场可以等效为入射场u0和散射场u1的叠加；
u＝u0+u1 (2)
其中入射场满足均匀介质方程：
忽略高阶无穷小量
联合式(1)、(2、(3)可以得到散射场与入射场的关系方程：
进一步利用格林函数和叠加原理推导出非均匀区域内的散射声场为：
其中x′为非均匀介质区域V中的一点，G为入射波场u0对应的格林函数，R＝|x-x′|表示两点之间的距离；将非均匀区域离散化为更小的体积单元，进而通过时域积分方程法叠加所有子单元的声场得到非均匀介质的散射声场。
优选的，所述步骤1-3中利用偏移成像算法得到井周地质结构的标签成像图的具体操作为：
步骤1-3-1:通过步骤1-2中得到的波形数据，使用Kirchhoff方程将波形数据反向传播至偏移域，偏移域中某点ζ0在t0时刻的反向传播声场表达式为：
式中，
表示偏移域Ω的边界，其中
为Ω外法线方向的导数；
步骤1-3-2:用高斯束GGB(ζ,t；ζ0,t0)叠加的渐近形式来替换式(6)中格林函数，得到偏移域中某点ζ0在t0时刻反向传播的波场表达式：
井中ζs位置处声源激发的直达波场的波动方程为：
式中，f(t)为测井仪器的声源激发函数，一般为雷克子波；通过与上述反向传播波场相同的推导步骤得到直达波场的表达式为：
式中，Re表示复数取实部，fF(ω)为激发函数的频域表示，GGB(ζ,ζs；ω)表示从声源点ζs到偏移域内点ζ的高斯束簇，为波动方程中格林函数的渐近解；
格林函数的渐近解GGB(ζ,ζ0；ω)是由一系列由点ζ0出射，并且经过成像点ζ的高斯束组成，其中单条高斯束射线在射线中心坐标系下的表达式为：
s为射线的弧长，n为射线外一点到射线的距离，υ(s)为射线上的速度，p,q为动力学参数，τ(s)为中心射线上的到时；通过对不同出射方向的高斯束进行叠加，得到格林函数的渐近解：
式中，Φ(θ；ω)是出射角为θ的高斯束的初始幅度，公式如下所示：
将上述频域表示的格林函数变换到时域：
步骤1-3-3:利用式(10)，对高斯束进行求解后，再利用式(11)和式(9)即可计算直达波场，进而利用式(13)和式(7)计算反向传播波场，接下来计算两种波场的互相关进行成像，公式如下：
I(ζ；ζs)＝∫dt0uD(ζ,t0；ζs)u(ζ,t0) (14)
上式通过计算直达波场uD(ζ,t0；ζs)与反向传播波场u(ζ,t0)之间的相关值来对井外地质结构进行成像。
本发明的有益效果：1、通过基于图像处理的角度对声波远探测成像图进行直接降噪，构建降噪神经网络模型，训练端到端的降噪神经网络，不需要复杂的参数设置，在保留成像图边缘信息的同时，自动快速地移除成像图中干扰噪声，提高了成像图的保存质量。
2、通过Born近似的散射声场计算方法和偏移成像方法建立仅包含反射体的标签成像图数据集，能够支持深度卷积神经网络的有效训练，避免了繁琐费时的现场数据搜集及标注过程。
附图说明
图1为本发明的方法流程图。
图2为建立的裂缝型地层模型。
图3为地层模型的偏移成像图。
图4为地层模型的偏移成像图中添加了白高斯噪声后的噪声成像图。
图5为添加了白高斯噪声后的噪声成像图通过本发明方法处理后得到的降噪后成像图。
图6为构建的降噪网络模型示意图。
图7为验证集上损失函数的迭代优化曲线图。
图8为验证集上信噪比的迭代优化曲线图。
图9为训练后网络用于现场成像图的处理结果图。
具体实施方式
以下结合附图对本发明做进一步概况。
如图1所示，一种声波远探测成像图的降噪方法，包括以下步骤：
步骤1：准备数据集；在均匀地层的基础上叠加非均匀介质构建地质模型，通过计算井外非均匀介质的散射效应评价远探测接收波形，利用偏移成像算法处理时域波形得到标签成像图，随机添加白高斯噪声得到噪声成像图，随机生成并处理大量地质模型，建立大规模模拟数据集。
步骤2：构建降噪网络模型；基于卷积神经网络建立全分辨率的深度网络模型，设定其中的模型深度、卷积核尺寸、步长和维度、激活函数类型、批归一化层。
步骤3：数据预处理；将步骤一中建立的数据集分为训练、验证和测试集，将成像图进行归一化，并采用随机偏移、裁剪、缩放、翻转等方式进行数据增广。
步骤4：训练网络模型；设定批次采样数量，在训练集中进行随机采样，将噪声成像图和对应的噪声等级值输入给步骤二中构建的网络模型，联合网络的输出与对应的标签成像图计算L1损失函数，通过随机梯度下降法对网络进行优化，直至达到预定的迭代次数，或者验证集上的评价指标变化达到预设阈值。
步骤5：现场成像图处理；将步骤4中训练结束的模型参数进行固化，并根据模拟数据集中成像图的尺寸对现场成像图进行滑动重叠切片，将切片输入给网络进行降噪，最终反向合并各切片对应输出，得到最终的降噪后成像图。
结合实施例进一步说明：
步骤1：准备数据集：
步骤1-1：设定均匀地层为背景介质，尺寸为20×20m,纵波速度为6220m/s,横波速度为3208m/s，密度为2700kg/m3；裂缝型非均匀介质的纵波速度为4010m/s,横波速度为2102m/s，密度为2650kg/m3；将非均匀区域叠加在均匀地层之上，构建地质模型；构建地质模型如图2所示；
步骤1-2：利用Born近似对局部非均匀介质引起的平面波散射场进行处理，首先通过求解均匀介质的波动方程评价入射波场，然后对非均匀区域进行离散化，对各网格声场响应进行叠加评价散射波场，并通过时域积分方程法评价所有时刻的波形数据；
步骤1-2-1：当井外地层中存在传播声速变化的非均匀介质时，井中声源激发的弹性波在该区域内将会改变其原始传播方向，产生散射声场；
步骤1-2-2：在弱散射条件下，非均匀介质速度变化引起的局部不均匀性，散射波的传播方程由入射波场与介质的非均匀性确定：
其中u为标量位移势，V0为非均匀介质的背景速度；
步骤步骤1-2-3：基于Born近似原理，总声场可以等效为入射场u0和散射场u1的叠加；
u＝u0+u1 (7)
其中入射场满足均匀介质方程：
忽略高阶无穷小量
联合式(1)、(2)、(3)可以得到散射场与入射场的关系方程：
进一步利用格林函数和叠加原理推导出非均匀区域内的散射声场为：
其中x′为非均匀介质区域V中的一点，G为入射波场u0对应的格林函数，R＝|x-x′|表示两点之间的距离；将非均匀区域离散化为更小的体积单元，进而通过时域积分方程法叠加所有子单元的声场得到非均匀介质的散射声场。
步骤1-3：根据均匀地层声学参数构建速度模型，利用偏移成像算法评价得到井周地质结构的标签成像图，如图3所示，其中反射体与地质模型基本相互对应；
所述步骤1-3中利用偏移成像算法得到井周地质结构的标签成像图的具体操作为：
步骤1-3-1:通过步骤1-2中得到的波形数据，使用Kirchhoff方程将波形数据反向传播至偏移域，偏移域中某点ζ0在t0时刻的反向传播声场表达式为：
式中，
表示偏移域Ω的边界，其中
为Ω外法线方向的导数；
步骤1-3-2:用高斯束GGB(ζ,t；ζ0,t0)叠加的渐近形式来替换式(6)中格林函数，得到偏移域中某点ζ0在t0时刻反向传播的波场表达式：
井中ζs位置处声源激发的直达波场的波动方程为：
式中，f(t)为测井仪器的声源激发函数，一般为雷克子波；通过与上述反向传播波场相同的推导步骤得到直达波场的表达式为：
式中，Re表示复数取实部，fF(ω)为激发函数的频域表示，GGB(ζ,ζs；ω)表示从声源点ζs到偏移域内点ζ的高斯束簇，为波动方程中格林函数的渐近解；
格林函数的渐近解GGB(ζ,ζ0；ω)是由一系列由点ζ0出射，并且经过成像点ζ的高斯束组成，其中单条高斯束射线在射线中心坐标系下的表达式为：
s为射线的弧长，n为射线外一点到射线的距离，υ(s)为射线上的速度，p,q为动力学参数，τ(s)为中心射线上的到时。通过对不同出射方向的高斯束进行叠加，得到格林函数的渐近解：
式中，Φ(θ；ω)是出射角为θ的高斯束的初始幅度，公式如下所示：
将上述频域表示的格林函数变换到时域：
步骤1-3-3:利用式(10)，对高斯束进行求解后，再利用式(11)和式(9)即可计算直达波场，进而利用式(13)和式(7)计算反向传播波场，接下来计算两种波场的互相关进行成像，公式如下：
I(ζ；ζs)＝∫dt0uD(ζ,t0；ζs)u(ζ,t0) (14)
上式通过计算直达波场uD(ζ,t0；ζs)与反向传播波场u(ζ,t0)之间的相关值来对井外地质结构进行成像。
步骤1-4：给标签成像图添加等级为35的全局白高斯噪声，得到噪声成像图，如图4所示，反射体基本被严重的噪声所淹没；
步骤1-5：重复步骤1-1、步骤1-2、步骤1-3和步骤1-4，得到1万组成对的标签和噪声成像图，建立大规模的理论模拟数据集。
步骤2：构建如图6所示的降噪网络模型：
步骤2-1：通过间隔取值的方式的对输入噪声成像图进行可逆降采样，生成四组尺寸降低一倍的含噪声子图；
步骤2-2：建立与含噪声子图相同尺寸的二维噪声等级图，并填充所添加的白高斯噪声的等级值，将含噪声子图与噪声等级图进行拼接；
步骤2-3：利用尺寸为3×3、步长为1×1、维度为5的卷积核对网络输入进行处理，并采用ReLU激活函数进行非线性变换，然后通过尺寸为3×3、步长为1×1、维度为64的卷积核、ReLU激活函数和批归一化层构建卷积层，并连续累加14组卷积层在同一分辨率进行特征提取，充分保存图像中的位置信息，然后利用尺寸为3×3、步长为1×1、维度为4的卷积核将提取特征转化为降噪后的子图，通过可逆降采样的反操作合并各子图得到降噪后成像图。
步骤3：数据预处理：
步骤3-1：按照75％、15％、10％的比例将数据集分解为训练集、验证集和测试集；
步骤3-2：对图像进行归一化：
其中μ表示数据的均值，σ为数据的标准差；
步骤3-3：随机选择偏移、裁剪、缩放、翻转等方式对数据进行增广。
步骤4：训练神经网络：
步骤4-1：按照批次采样数为64从训练集中进行随机采样，将噪声成像图和噪声等级值输入给图3中构建的降噪神经网络；
步骤4-2：联合网络输出与对应的标签数据计算L1损失函数：
其中N为批次采样值，p为网络输出，y为标签成像图，|·|为绝对值运算符。
步骤4-3：利用随机梯度下降法将计算的L1损失函数反向传播，更新神经网络参数；
步骤4-4：重复步骤4-1至步骤4-3对网络进行连续优化，直至训练集中所有数据均被采样完成，对验证集数据进行处理，计算L1损失函数与信噪比指数：
SNR＝10log10(‖p‖2/‖p-y‖2)
步骤4-5：重复步骤4-4，直至迭代次数达到200，或者5次迭代过程中验证集上的信噪比指数上升小于0.01。如图7所示：纵轴表示L1损失函数，横轴表示迭代次数，L1损失函数从初始的0.015逐步降低至0.0039，如图8所示：纵轴表示信噪比指标，横轴表示迭代次数，给出了对应的信噪比指标，经过78次迭代过程后，信噪比提升至38.4，表明网络得到了良好优化。
步骤4-6：将训练好的网络应用于测试数据集，评价网络的泛化性能，如图5中所示，经过网络处理后，图中的噪声基本被完全移除，图中的反射体的边缘得到良好保存，证明了网络的降噪效果。
步骤5：现场成像图降噪：
步骤5-1：将步骤4中训练结束的模型参数进行固化和保存；
步骤5-2：对如图9中第二道所示的现场成像图进行处理，首先按照256×256的尺寸对成像图进行滑动切片，并且保持64×64的重叠区域；
步骤5-3：设定噪声等级值为30，与分别归一化后的切片一起输入给降噪网络，将各输出进行反归一化后重新拼接，对于重叠区域则采用取平均值的方式，最终得到降噪后的成像图，如图9中第三道所示，其中的背景噪声基本被移除，反射体变得清晰可见，视觉效果得到提升。
工作原理：声波远探测测井技术利用井中声源向井外地层激发弹性波，经井外非均质体反射后被井中接收器采集，进而通过偏移成像算法得到井外地质结构的高分辨率成像图，为了移除成像图中干扰噪声，首先联合基于Born近似的散射声场计算方法和偏移成像方法建立仅包含反射体的标签成像图数据集，然后通过添加不同标准差的高斯噪声模拟实际成像图中的背景噪声，针对井外地质结构尺度、位置的随机分布特性，利用神经网络充分提取成像图数据集中反射体的特征信息，在此基础上智能移除背景噪声，良好保留反射体的边缘信息，为后续的储层识别提供更清晰准确的井周成像图。
以上所述仅为本发明的实施方式而已，并不用于限制本发明。对于本领域技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原理内所作的任何修改、等同替换、改进等，均应包括在本发明的权利要求范围之内。
Claims (8)
1.一种声波远探测成像图的降噪方法，其特征在于：所述该方法包括以下步骤：
步骤1：准备数据集；
步骤2：构建降噪网络模型；
步骤3：数据预处理；
步骤4：训练网络模型；
步骤5：现场成像图处理。
2.根据权利要求1所述的方法，其特征在于：所述步骤1中准备数据集的具体操作为：
步骤1-1：利用均匀地层作为背景介质，设定其尺寸、纵横波速度及密度值，设定非均匀介质的纵横波速度和密度值，并利用随机函数确定其尺寸、形状、位置，将非均匀区域叠加在均匀地层之上，构建地质模型；
步骤1-2：通过求解均匀介质的波动方程评价入射波场，然后对非均匀区域进行离散化，对各网格声场响应进行叠加评价散射波场，通过时域积分方程法评价所有时刻的接收波形数据；
步骤1-3：利用偏移成像算法将中得到的时域波形向倾斜同相轴下倾方向归位，得到井周地质结构的标签成像图；
步骤1-4：给标签成像图随机添加不同等级的白高斯噪声，得到噪声成像图；
步骤1-5：重复步骤1-1、1-2、1-3、1-4，得到大量成对的标签和噪声成像图，建立大规模的理论模拟数据集。
3.根据权利要求1所述的方法，其特征在于：所述步骤2中构建降噪网络模型采用深层卷积神经网络，构建深层卷积神经网络包括以下步骤：
步骤2-1：通过间隔取值的方式的对尺寸为M×N×1的输入噪声成像图进行可逆降采样，生成尺寸为
的含噪声子图，其中最后一维表示通道；
步骤2-2：建立与含噪声子图相同尺寸
的二维噪声等级图，并全部填充为所添加的白高斯噪声的等级值，将含噪声子图与噪声等级图进行拼接，得到尺寸为
的网络输入；
步骤2-3：利用尺寸为3×3、步长为1×1、维度为5的卷积核对网络输入进行处理，并采用ReLU激活函数进行非线性变换，然后通过尺寸为3×3、步长为1×1、维度为64的卷积核、ReLU激活函数和批归一化层构建卷积层，并连续累加14组卷积层在同一分辨率进行特征提取，充分保存图像中的位置信息；利用尺寸为3×3、步长为1×1、维度为4的卷积核将提取特征转化为降噪后的子图，通过可逆降采样的反操作合并各子图得到降噪后成像图。
4.根据权利要求1所述的方法，其特征在于：所述步骤3中数据预处理的具体操作为：
步骤3-1：按照75％、15％、10％的比例将数据集分解为训练集、验证集和测试集；
步骤3-2：对图像进行归一化：
其中μ表示数据的均值，σ为数据的标准差；
步骤3-3：随机选择偏移、裁剪、缩放、翻转的方式对数据进行增广。
5.根据权利要求1所述的方法，其特征在于：所述步骤4中训练网络模型的具体操作为：
步骤4-1：按照设定的批次采样数从训练集中进行随机采样，将噪声成像图和噪声等级值输入给步骤2中构建的降噪神经网络模型；
步骤4-2：联合网络输出与对应的标签数据计算L1损失函数：
其中N为批次采样值，p为网络输出，y为标签成像图，|·|为绝对值运算符；
步骤4-3：利用随机梯度下降法将计算的L1损失函数反向传播，更新神经网络参数；
步骤4-4：重复步骤4-1至步骤4-3对网络进行连续优化，直至训练集中所有数据均被采样完成，对验证集数据进行处理，计算L1损失函数与信噪比指数：
SNR＝10log10(‖p‖2/‖p-y‖2)；
其中，SNR为信噪比指数；
步骤4-5：重复步骤4-4，直至迭代次数达到200，或者5次迭代过程中验证集上的信噪比指数上升小于0.01；
步骤4-6：将训练好的网络应用于测试数据集，评价网络的泛化性能。
6.根据权利要求1所述的方法，其特征在于：所述步骤5中现场成像图处理的具体操作为：
步骤5-1：将步骤4中训练结束的模型参数进行固化和保存；
步骤5-2：首先按照理论模拟成像图的尺寸对成像图进行滑动切片，在边界处存在重叠区域；
步骤5-3：将分别归一化后的切片与设定的噪声等级值输入给降噪网络，将各输出进行反归一化后重新拼接，对于重叠区域则采用取平均值的方式，最终得到降噪后的成像图。
7.根据权利要求2所述的方法，其特征在于：所述步骤1-2中通过求解均匀介质的波动方程评价入射波场的具体操作为：
步骤1-2-1：当井外地层中存在传播声速变化的非均匀介质时，井中声源激发的弹性波在该区域内将会改变其原始传播方向，产生散射声场；
步骤1-2-2：在弱散射条件下，非均匀介质速度变化引起的局部不均匀性，散射波的传播方程由入射波场与介质的非均匀性确定：
其中u为标量位移势，V0为非均匀介质的背景速度；
步骤1-2-3：基于Born近似原理，总声场可以等效为入射场u0和散射场u1的叠加；
u＝u0+u1(2)
其中入射场满足均匀介质方程：
忽略高阶无穷小量
联合式(1)、(2)、(3)可以得到散射场与入射场的关系方程：
进一步利用格林函数和叠加原理推导出非均匀区域内的散射声场为：
其中x′为非均匀介质区域V中的一点，G为入射波场u0对应的格林函数，R＝|x-x′|表示两点之间的距离；将非均匀区域离散化为更小的体积单元，进而通过时域积分方程法叠加所有子单元的声场得到非均匀介质的散射声场。
8.根据权利要求2所述的方法，其特征在于：所述步骤1-3中利用偏移成像算法得到井周地质结构的标签成像图的具体操作为：
步骤1-3-1:通过步骤1-2中得到的波形数据，使用Kirchhoff方程将波形数据反向传播至偏移域，偏移域中某点ζ0在t0时刻的反向传播声场表达式为：
式中，
表示偏移域Ω的边界，其中
为Ω外法线方向的导数；
步骤1-3-2:用高斯束GGB(ζ,t；ζ0,t0)叠加的渐近形式来替换式(6)中格林函数，得到偏移域中某点ζ0在t0时刻反向传播的波场表达式：
井中ζs位置处声源激发的直达波场的波动方程为：
式中，f(t)为测井仪器的声源激发函数，一般为雷克子波；通过与上述反向传播波场相同的推导步骤得到直达波场的表达式为：
式中，Re表示复数取实部，fF(ω)为激发函数的频域表示，GGB(ζ,ζs；ω)表示从声源点ζs到偏移域内点ζ的高斯束簇，为波动方程中格林函数的渐近解；
格林函数的渐近解GGB(ζ,ζ0；ω)是由一系列由点ζ0出射，并且经过成像点ζ的高斯束组成，其中单条高斯束射线在射线中心坐标系下的表达式为：
s为射线的弧长，n为射线外一点到射线的距离，υ(s)为射线上的速度，p,q为动力学参数，τ(s)为中心射线上的到时；通过对不同出射方向的高斯束进行叠加，得到格林函数的渐近解：
式中，Φ(θ；ω)是出射角为θ的高斯束的初始幅度，公式如下所示：
将上述频域表示的格林函数变换到时域：
步骤1-3-3:利用式(10)，对高斯束进行求解后，再利用式(11)和式(9)即可计算直达波场，进而利用式(13)和式(7)计算反向传播波场，接下来计算两种波场的互相关进行成像，公式如下：
I(ζ；ζs)＝∫dt0uD(ζ,t0；ζs)u(ζ,t0) (14)
上式通过计算直达波场uD(ζ,t0；ζs)与反向传播波场u(ζ,t0)之间的相关值来对井外地质结构进行成像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210845762.6A CN115170428A (zh) | 2022-07-18 | 2022-07-18 | 一种声波远探测成像图的降噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210845762.6A CN115170428A (zh) | 2022-07-18 | 2022-07-18 | 一种声波远探测成像图的降噪方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115170428A true CN115170428A (zh) | 2022-10-11 |
Family
ID=83495567
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210845762.6A Pending CN115170428A (zh) | 2022-07-18 | 2022-07-18 | 一种声波远探测成像图的降噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115170428A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115576014A (zh) * | 2022-10-26 | 2023-01-06 | 江苏科技大学 | 一种基于声波远探测成像的裂缝型储层智能识别方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108765319A (zh) * | 2018-05-09 | 2018-11-06 | 大连理工大学 | 一种基于生成对抗网络的图像去噪方法 |
| CN111860273A (zh) * | 2020-07-14 | 2020-10-30 | 吉林大学 | 基于卷积神经网络的磁共振地下水探测噪声抑制方法 |
| CN113156515A (zh) * | 2021-04-16 | 2021-07-23 | 中国石油大学(华东) | 一种声波远探测成像降噪处理方法和装置 |
-
2022
- 2022-07-18 CN CN202210845762.6A patent/CN115170428A/zh active Pending
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108765319A (zh) * | 2018-05-09 | 2018-11-06 | 大连理工大学 | 一种基于生成对抗网络的图像去噪方法 |
| CN111860273A (zh) * | 2020-07-14 | 2020-10-30 | 吉林大学 | 基于卷积神经网络的磁共振地下水探测噪声抑制方法 |
| CN113156515A (zh) * | 2021-04-16 | 2021-07-23 | 中国石油大学(华东) | 一种声波远探测成像降噪处理方法和装置 |
Non-Patent Citations (3)
| Title |
|---|
| KAI ZHANG等: "FFDNet: Toward a Fast and Flexible Solution for CNN based Image Denoising", IEEE TRANSACTIONS ON IMAGE PROCESSING, vol. 27, no. 9, pages 4610 - 4611 * |
| 李盛清: "声波远探测成像处理方法及地质应用研究", 中国博士学位论文全文数据库基础科学辑, pages 2 * |
| 黄建平 等: "格林函数高斯束逆时偏移", 石油地球物理勘探, vol. 49, no. 01, pages 102 - 103 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115576014A (zh) * | 2022-10-26 | 2023-01-06 | 江苏科技大学 | 一种基于声波远探测成像的裂缝型储层智能识别方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Van den Ende et al. | A self-supervised deep learning approach for blind denoising and waveform coherence enhancement in distributed acoustic sensing data | |
| CN106772583A (zh) | 一种地震绕射波分离方法和装置 | |
| CN102298156A (zh) | 用于反虚反射地震数据的方法和装置 | |
| CN103675910B (zh) | 一种水陆检波器地震数据标定因子反演方法 | |
| WO2009002001A1 (en) | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging | |
| MX2014003060A (es) | Sistemas y metodos de filtracion de dominio de frecuencia y discriminacion de dominio de espacio-tiempo de datos sismicos. | |
| US20110199858A1 (en) | Estimating internal multiples in seismic data | |
| CN113015926A (zh) | 无源地震成像 | |
| RU2570827C2 (ru) | Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников | |
| EP3211594B1 (en) | Seismic modeling system providing seismic survey data inpainting based upon suspect region boundary comparisons and related methods | |
| EP3217354B1 (en) | Seismic modeling system providing seismic survey data frequency domain inpainting and related methods | |
| CN116520419A (zh) | 一种热流体裂缝通道识别方法 | |
| CN104199088B (zh) | 一种提取入射角道集的方法及系统 | |
| CN115170428A (zh) | 一种声波远探测成像图的降噪方法 | |
| US11143769B2 (en) | Seismic modeling system providing seismic survey data spatial domain exemplar inpainting and related methods | |
| AU2016220145B2 (en) | Black hole boundary conditions | |
| KR101318994B1 (ko) | 복수의 가중치를 이용한 지하매질 구조 추정 방법 및 장치 | |
| CN115576014A (zh) | 一种基于声波远探测成像的裂缝型储层智能识别方法 | |
| US10338253B2 (en) | Method of suppressing spectral artefacts of wavefield decomposition caused by imperfect extrapolation | |
| CN114740528A (zh) | 一种超微分拉普拉斯块约束的叠前多波联合反演方法 | |
| US10338250B2 (en) | Method of removing incoherent noise | |
| CN116088049B (zh) | 基于子波变换的最小二乘逆时偏移地震成像方法及装置 | |
| CN107664769B (zh) | 一种角度域共成像点道集提取方法和装置 | |
| CN116011306A (zh) | 基于平面震源地震数据的速度建模方法、装置及电子设备 | |
| CN117665930A (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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20221011 |
|
| RJ01 | Rejection of invention patent application after publication |