各向异性扩散滤波远探测声波测井图像降噪方法*

傅艳莉 李 超 陈 浩 曹雪砷

(1 中国科学院大学 北京 100049)

(2 中国科学院声学研究所 声场声信息国家重点实验室 北京 100190)

(3 北京市海洋深部钻探测量工程技术研究中心 北京 100190)

0 引言

近年来，随着石油、天然气及矿藏勘探和开发需求的增长及测井技术的发展，声波远探测测井成为井旁构造体识别的热点技术[1]。常规测井项目探测井外信息较浅，难以反映距井较远距离的地质异常体的发育情况。声波远探测技术可以对井周围数十米甚至上百米范围进行探测，利用地层反射信息进行井周成像，探测到地层中裂缝、断层等，大大丰富了研究领域，已经发展成一种深部复杂油气勘探领域中不可缺少的先进技术[2-3]。然而由于反射波幅度低、受幅度强的井筒直达波干扰等原因，反射波图像具有低信噪比的特点，噪声主要是残留的井孔模式波和随机噪声(包括地层噪声和滤波处理过程中引入的噪声)，这些噪声大都是非相干的。传统的图像去噪方法，如均值滤波、中值滤波等，在去除噪声的同时会造成裂缝、断层等地质体结构的模糊。因此，为了解决去噪同时兼顾构造边缘及细节的问题，前人引入了扩散滤波处理远探测测井图像[4]，取得了较好效果。但目前尚缺少对该方法适用性及其特点的系统性研究和分析，而采用各向异性扩散滤波对地层数据进行处理，在去除噪声的同时可以更好地保留裂缝、断层等地质体结构[5]。

基于扩散方程的噪声衰减方法是20 世纪90年代提出的一种数字图像处理方法，是图像降噪的研究热点。与传统的线性滤波方法相比，扩散滤波在去噪的同时保留图像的边缘细节，广泛应用于图像降噪领域的研究。其中最经典的是Perona和Malik提出的各向异性扩散模型[6]。但是P-M模型的边缘保持效果不明显，后来又出现了许多改进模型，其中Weickert[7]模型用结构张量的形式表示扩散系数，使得各向异性扩散滤波广泛应用于地震勘探的图像处理中。Bakker[8]提出将基于图像的各向异性扩散滤波应用到地震资料的预处理，并且提出几何张量结构算法，用地层的几何结构张量中包含的丰富地层信息，来描述地层的内部结构，计算检测断裂系统。Fehmers等[9]把各向异性扩散滤波引入地震信息资料结构信息的边缘保持。Lavialle等[10]提出地震断层保持的扩散滤波方法，使得处理之后的地震数据中断层信息更加突出。孙夕平等[11]、王绪松等[12]将二阶导数引入扩散张量，探索了一致增强扩散滤波技术在二维地震剖面的保边滤波应用。杨培杰等[13]提出一种方向性边界增强技术，对于断层信息起到了很好的保持作用。张尔华等[14]对三维地震资料非线性各向异性扩散滤波方法进行了初步探讨，通过控制在断层等不连续区域的滤波程度，使得各向异性扩散滤波器具有较好的保边处理性能。各向异性扩散滤波在地震图像的处理中已经有了广泛的应用，本文研究该降噪方法在声波远探测处理的适用性。

本文首先分析了远探测声波测井的残余噪声问题，通过梳理近年来各向异性扩散滤波的发展脉络，从基本原理入手并针对其在远探测图像处理中的应用研究展开分析，对比具有不同扩散张量的滤波方法及滤波参数的处理效果，建立了适合于远探测声波测井信号的处理方法流程，最后通过实际数据的处理分析验证了方法的效果，并对该技术今后改进和发展方向进行了展望与讨论。

1 远探测声波测井残余噪声分析

由于反射波幅度低、受幅度强的井筒直达波干扰等原因，反射波图像具有低信噪比的特点，原始波形经过反射波分离处理后，广义上的噪声主要是残留的井孔模式波和随机噪声(包括地层噪声和滤波处理过程中引入的噪声)，除了残余模式波外，噪声大都是非相干的，对后续的处理和地质体识别造成干扰。以深部地层采集的偶极远探测实际数据为例，如图1(a)所示，原始的反射波分离后的剖面信噪比较低，部分反射体特征淹没在较强的噪声中，除了10 ms以内的在轴向上表现为非相干的残余弯曲波外，图中可见两种斜率的反射波，如黑色圆圈和红色圆圈所示，其特征表现为在时间上延伸较长，顶部和井孔模式波重合。为了说明噪声特征并区别其来源是否为井外反射，取其中3 个深度点阵列波形进行相关速度分析。图1(b)为对图像6350 m深度的8道波形(见图1(a)区域A)做的相干分析，由于未见明显反射，分析表明其是非相干的随机噪声。图1(c)为对深度6245 m 的阵列波形(见图1(a)区域B)进行相关分析，正向阵列(从1 到8)未见相干性，但反向阵列(从8 到1)得到其视速度在1400 m/s 左右，由于反射波在接收阵列上的入射角较小，一般具有较高的视速度，因此确定其为沿井轴反向传播的反射斯通利波(下行)。同理，图1(d)是对另一斜率的6310 m 深度波形(见图1(a)中区域C)相干分析的结果，得到其下行视速度为3100 m/s 左右，验证其为沿井轴传播的反射横波。

图1 反射波分离后剖面及3 个深度点的相干分析Fig.1 Reflection wave separation profile and coherence analysis of three depth points

以上分析表明，随机噪声、残余模式波和井孔模式波反射存在于远探测声波测井的波形中，尤其是沿井轴传播的反射干扰常被误以为是井外反射。这些噪声的存在，使获取的图像难以表示真正的井外反射体特征。而一般而言，经过叠加后，反射信号得到一定程度的增强，这些噪声特征会减弱，总体表现为非相干性，但其仍然会掩盖部分真实反射波。因此采取适当的滤波方法进一步压制噪声并增强反射信号是实际中需要解决的关键难题。针对该问题，本文引入各向异性扩散滤波方法来进行处理分析。

2 各向异性扩散滤波

2.1 各向异性扩散滤波原理

各向异性扩散滤波是一种基于偏微分方程的图像处理方法。在物理学中，扩散定义是指物质分子从高热量区域向低热量区域转移，最终达到平衡状态。当把扩散过程应用到图像处理中，研究对象就变成了图像的像素值。最简单的扩散为线性各向同性扩散，Koenderink[15]证明图像与高斯核函数的卷积与热传导方程是一致的，将物理学的热扩散偏微分方程应用到图像处理中可以表示为

式(1)中，U是待处理的数据；a为扩散系数，在各向同性扩散中a是一个常数；t为扩散时间。由于上述滤波器的扩散系数为常量，即对图像各个方向的扩散强度是一样的，不能起到保护图像边缘和细节信息的作用。为了克服这个缺点，Perona等[6]提出了非线性的各向异性扩散滤波模型：

式(2)中，将扩散系数设定为关于图像梯度的非负单减函数g(‖∇U ‖)，满足g(0)=1，g(∞)=0。但是P-M 模型由于选取的扩散系数是一个标量，没有考虑到图像的纹理特征，使得扩散结果存在阶梯效应。

Weickert[7]将结构分析引入扩散滤波中，将扩散系数变为结构张量，使扩散随方向变化而变化，在去噪的同时保护了图像的线性纹理特征。构建了如下扩散模型：

式(4)中，*为卷积运算符，G为高斯核函数，σ为噪声尺度，用以控制滤波过程中需要保护对象的最小尺寸，ρ为整合尺度，一般整合尺度应大于噪声尺度。结构张量Sρ是2×2 的对称半正定矩阵，对其进行特征分解：

式(5)中，v1和v2为梯度结构张量的两个特征向量,v1是梯度变化最大的方向，v2是垂直于梯度变化的方向；λ1和λ2为梯度结构张量相对应的非负特征值，且λ1≥λ2，特征值反映了信号在局部特征方向的变化强度。

扩散的方向由结构张量的特征向量决定，要保证扩散张量的特征向量与结构张量的特征向量相一致：

式(6)中，μ1和μ2是两个扩散方向的扩散强度，主要有3 种设置方法：相干增强扩散模型[16]、边缘增强扩散模型[17]和混合扩散模型[18]。

相干增强扩散模型在执行过程中，会使得待处图像的每个像素都沿着相干方向进行扩散，它的特征值取值与结构张量Sρ的两个特征值λ1和λ2的关系有关，设定为

式(7)中，参数α为一很小的正实数，控制扩散滤波过程中对不连续结构特征的保持程度，一般取0.001，使图像梯度变化较大方向的扩散强度很小，保护图像的结构信息；k表征图像不同方向的相干性，表示为k=(λ1-λ2)2；c1为阈值，通常取为1。

在图像相干性较高的区域，μ2约等于1，也就是说沿着特征向量v1方向的扩散强度很小，近似于不进行扩散，主要沿着特征向量v2的方向进行扩散；而在各向同性区域，v1和v2两个方向的扩散强度都很小，因此可以有效地保护图像的结构信息。

边缘增强扩散体现在沿边缘方向扩散大而垂直边缘方向扩散小，以保护图像边缘结构。其扩散张量的特征值设定为

式(8)中，c2和e为正实数，分别取3 和0.02；gradA为图像U两个方向梯度的平方和。μ2=1即在垂直于梯度方向扩散强度始终为1；而在梯度方向的扩散强度取决于图像两个方向梯度的平方和gradA,这个值越大，μ1越小，说明梯度方向的扩散越小。有效地识别图像边缘，对边缘起到保护与增强作用。

混合各向异性扩散是将相干增强扩散与边缘增强扩散结合起来。其扩散张量的特征值设定为

式(11)中，h为一正实数，取0.5；其他变量与式(7)、式(8)中取相同值。

式(3)可以用有限差分进行求解，主要有两种离散方法：半隐式离散形式[16]和显式离散形式[16]。利用差分代替微分，实现对图像U对扩散时间的求导，得到非线性扩散方程的半隐式离散形式：

式(13)中，Un+1和Un分别是图像在第n+1 和第n次迭代得到的扩散结果；Δt为迭代步长，取值越小表示扩散越慢，为了保证图像计算精度和稳定性，一般取一个很小的正实数。

使用正向差分近似对式(3)离散化，这样可以直接从前一个时间级别计算出新的时间级别上的值，而不用求解线性或非线性方程，得到显式离散形式：

式(14)中，An*Un是div(D∇u)的离散形式，即将图像与一个时变与空变的模板An进行卷积。An如表1[17]所示，其中，a、b和c是扩散张量中的3 个时、空变的值。

表1 卷积模板Table 1 Convolutional template

关于离散化方式，在同样的精度要求下，显式有限差分方法比半隐式格式更加有效；而半隐式扩散可以采用较大的步长，稳定性要求比显式低。

2.2 三种扩散模型的对比分析

相干增强扩散模型将μ1设置为一个极小的常数α，即在特征向量v1方向的扩散很小，而在v2方向的扩散取决于图像的相干性，图像相干性强的区域(图像边缘)μ2接近于1，该模型适用于图像噪声比较强的情况，用于增强线状纹理；边缘增强扩散模型正好相反，将μ2设置为常数1，即在垂直梯度方向的扩散始终保持最强，而在v1方向的扩散取决于图像两个方向梯度的平方和，该模型适用于图像噪声较弱的情况，用于平滑噪声，同时增强边缘；混合扩散模型是将以上两种模型结合起来，根据不同的噪声强度和图像特征选择不同的扩散模型。为了更形象地说明3 种扩散张量模型的差异，以指纹图像为例对其进行进一步对比分析。如图2(c)、图2(f)中红色箭头所示位置，相比于其他两种扩散模型，相干增强扩散滤波后间断图像的连续性增强，图像的纹理特征得到加强。

图2 指纹图像的滤波效果Fig.2 The filtering effect of fingerprint image

2.3 算法实现及参数选择

根据以上分析，设计的图像扩散迭代算法为

(1)输入原始地层图像U0，作为迭代算法的初始值，选择迭代次数K，设定当前迭代次数k=0；

(2)对图像做高斯平滑处理，得到Uσ；

(3)求Uσ的梯度结构张量，并对其做高斯平滑处理，得到Sρ；

(4)对Sρ做特征分解，得到两个特征值λ1和λ2及其对应的特征向量v1和v2；

(5)选择一种扩散张量模型计算扩散张量D；

(6)选择一种迭代方式得到迭代k+1次结果；

(7)如果k ＜K，则令k=k+1，返回步骤(2)；否则，结束迭代；

(8)输出滤波结果。

在步骤(2)中对高斯图像做尺度为σ的高斯卷积，来控制滤波过程中要保护对象的最小尺寸；在步骤(3)中对梯度结构张量做尺度为ρ的高斯卷积，可以消除噪声对计算结果的干扰。两个尺度参数的选择要根据待处理的数据特征确定，但要保证ρ ＞σ，当地层数据的裂缝结构比较丰富时，可以将两个尺度参数取小一点；当数据受到噪声干扰比较严重时，可以将尺度参数取大一点。在步骤(5)中利用梯度结构张量的特征值和特征向量构建扩散张量D，根据图像特性选择合适的扩散张量；在步骤(6)步中根据要处理的数据选择合适的迭代方法。迭代次数的选取要根据实际数据和预期的处理效果选定，保证处理后的数据信噪比有很大提高，裂缝得到很好的保持。同时，为了保证图像的收敛性，迭代步长Δt应该尽量小。

3 模拟数据处理分析

为了说明各向异性扩散滤波去除噪声同时保持图像细节的作用，利用模拟裂缝数据进行验证。原始模型数据为井外两列平行的裂缝的声波远探测偏移成像结果(图3(a))，在原始数据模型加入强噪声来模拟实际情况，噪声是使信噪比为0.1 的高斯噪声(图3(b))。可以看到，由于噪声的影响，裂缝特征变得模糊，给裂缝识别造成干扰。

图3(c)是用均值滤波方法处理的，在滤除噪声的同时，裂缝边缘变得不清晰，部分裂缝被消除，不利于裂缝结构的识别与解释。图4 是利用相干增强扩散张量模型、边缘增强扩散张量模型和混合扩散张量模型处理加噪数据得到的结果，采用的都是显式扩散。对比3 种模型的处理效果可以看出，3 种模型都能滤除图像中的大部分随机噪声，同时裂缝结构得到很好的显示。在深度15 m 附近的右边裂缝处，相干增强扩散张量模型比其余两种扩散张量模型更加清晰，处理效果更好，信噪比更高。为了更加清晰地看出3 幅处理后图像的区别，将其分别与原始图像做差得到图5(a)、图5(b)、图5(c)，可以看到5(a)右侧裂缝颜色更浅，说明裂缝结构得到了更好的保留，误差更小。而右侧裂缝相对左侧幅度较弱，说明相干增强扩散模型更适合于弱反射结构在强噪声环境下的处理，这是由于相对于梯度来说，相干性受噪声影响更小。图6(a)、图6(b)分别是采用相干增强扩散张量的半隐式扩散模型和显式扩散模型对加噪数据的处理结果，显式扩散模型的滤噪效果明显优于半隐式扩散模型。不仅有效滤除了图像中的噪声，保持了图像的分辨率，提高了信噪比，而且图像中的裂缝区域更加明显清晰，增强了裂缝特征的细节，具有很好的保持边缘的特性，增加了对裂缝的识别能力。一般而言，隐式差分比显示的稳定性要求要低，但是为了保持图像的细节特征，一般采用较小的时间步长，处理结果表明此时显式模型的效果要优于隐式的。

图3 模拟数据及均值滤波效果Fig.3 Simulation data and mean filtering result

图4 不同扩散张量模型的滤波效果Fig.4 Filtering results of different diffusion tensor models

图5 不同扩散张量模型的滤波图像与原图做差Fig.5 Filtering results of different diffusion tensor models are subtracted from the original image

图6 不同迭代模型的滤波效果Fig.6 Filtering results of different iteration models

为了定量地比较处理前后图像质量的改善效果，计算加噪图像与处理后图像的信噪比：

式(15)中，u0表示处理前的图像，u表示处理后的图像。选取不同的参数进行实验，并计算处理前后的信噪比和能量保持的关系。实际计算表明，处理后信噪比都增加了5～10 dB，能量保持60%以上,从图7 可以看出，在一定范围内信噪比随着迭代次数的增加而增加。相干增强显式扩散模型虽然信噪比增加较慢，但是随着迭代次数增加最终达到的信噪比最高。当迭代次数取35 时，图3 三种模型处理后的信噪比分别为9.03985 dB、8.6109 dB 和8.8968 dB。说明选取合适的参数，不仅可以抑制噪声还可以保留边界信息。

图7 不同扩散模型信噪比随迭代次数的变化Fig.7 The signal-to-noise ratio of different iteration models varies with the number of iterations

该实例说明，基于相干增强扩散张量的显式各向异性扩散滤波可以明显改善资料品质，有效识别裂缝结构。在处理中，选择σ= 2、ρ= 5、Δt= 0.1、K= 35，明显看出选择合适的参数与模型，各向异性扩散滤波处理之后的图像，信噪比得到提高，图像更加清晰，有利于裂缝结构的识别。因此，本文建立了相干增强扩散张量结合显式扩散模型的处理步骤，并依据该组合进一步对实际数据进行处理。

4 实际数据结果分析

下面是对远探测实测数据进行各向异性扩散滤波的一个实例。由于地下条件复杂，数据采集深度较大，在深部高温环境下换能器性能受到影响，信号质量下降，噪声增加。原始的反射波分离后的叠加剖面图像信噪比非常低,如图8(b)所示，虽然有部分反射体特征存在，但是均淹没在较强的噪声里面，与具有连续性的裂缝带相比，噪声在图像中的分布和大小不规则，更加杂乱且具备随机性和不连续性。图8(c)和图8(d)分别是采用相干增强扩散张量经过半隐式和显式各向异性扩散滤波之后的结果，3 幅图的色标范围一致。对比处理前后的剖面图像可以看出，处理之后的剖面相对于原始剖面，信噪比得到明显提高，裂缝带连续性也得到增强，特别是深度6310～6340 m 范围效果更明显。沿井壁附近传播并在时间上衰减较小的残留的反射斯通利波和反射横波(如图8(b)中黑色圆圈和红色圆圈所示)，由于经过反射波分离并叠加后其相干性较弱，也均得到了有效压制。图8(d)相比于图8(c)，更好地去除了反射斯通利波，如图8(b)第二个黑色圆圈所示。

图8 反射波分离结果剖面的各向异性扩散滤波效果对比Fig.8 Comparison of anisotropic diffusion filtering results of reflected wave separation result profile

图9 是另一深度段的反射波叠加剖面的处理结果，进一步验证了基于相干增强扩散张量的显式扩散模型在处理远探测实测数据上的优越性，信噪比更高，裂缝更加清晰，尤其是井壁附近的残留弯曲波噪声得到了很好的压制。本文中的数据降噪方法应用于实际数据，取得较好效果，改善了成像图的视觉质量，提高了成像中反射体的清晰度，增强了对井外构造的解释度。

图9 反射波分离结果剖面的各向异性扩散滤波效果对比Fig.9 Comparison of anisotropic diffusion filtering results of reflected wave separation result profile

5 结论

声波远探测作为一种有效的测井手段，在井外地质体探测和隐蔽油藏识别等方面有重要的作用。本文将图像处理中的各向异性扩散滤波技术引入声波远探测的图像处理中，经过对比将相干增强张量模型和显示扩散模型相结合，可以有效抑制深部地层中反射波图像中的非相干噪声，同时保持图像的地质结构信息，突出裂缝的边界特征。处理后的资料品质有较为明显的提高，对细微裂缝的分辨率更高，更有利于对小裂缝和断裂的识别，证实了各向异性扩散滤波在反射波图像处理中的有效性和实用性。各向异性扩散在强噪声干扰下的图像降噪方面可能会出现将有效信息误认为噪声的情况，这方面还有待进一步的完善和提高。为了进一步突显有效地质信息，可以将相干体技术和各向异性扩散滤波结合起来，进行属性提取分析，有望更好地表示井外反射体信息。