叠前非局部平均滤波压制随机噪音_胡新海.pdf

第 42 卷 第 5 期 煤田地质与勘探 Vol. 42 No.5 2014 年 10 月 COAL GEOLOGY 2. Institute of Porous Flow and Fluid Mechanics, Chinese Academy of Sciences, Langfang 065007, China; 3. Research Institute of Petroleum Exploration and Development-Langfang, Langfang 065007, China Abstract The nonlocal means has good denoising perance, but its application is newly developing in seismic data processing. The , using the structural redundancy of data, taking the small window with local structure and neighborhood as unit, conducts weighted arithmetic by using local structural similarity to enhance effective signals and to depress random noises. Aiming at huge amount of pre-stack seismic data, strong background noise and simple local structure, the original nonlocal means filters each point, conducts weighted calculation after calculating the weight coefficient of all points within data. Because of short points such as huge computation volume and poor adaptability to strong noise background, the original nonlocal means has been improved. Three modifications have been proposed for the nonlocal means algorithm. Firstly, the scan windows are divided with velocity spectrum; then, pre-selection of similar set is based on singular value decomposition in gradient domain; lastly, selection of self-adaptive filtering parameter is based on the scale of similar set. De-noising results for the test data demonstrate that the can effectively depress the random noise of seismic data. Key words pre-stack nonlocal means ; self-adaptive weighting; singular value decomposition in gradient domain; pre-selection; denoising 随机干扰严重降低了地震资料的信噪比，其频 谱很宽，无一定视速度，因而不能利用随机干扰和 有效波之间在频谱上的差异或传播方向上的差异 即视速度上的差异来进行压制[1]。 非局部平均滤波 Non-Local Means Filtering， NLMF是一种性能优异 的滤波方法，其利用图像具有的结构冗余度，以小 窗口或邻域为单元进行加权运算，增强有效信号的 结构压制随机噪音。相比常用的带通滤波、奇异值 分解、f-x 预测反褶积等方法，非局部滤波法是一种 全新的方法。原有的时间域和频率域滤波方法从本 质上都是对位置相邻的计算点加权平均的局部平均 滤波。非局部平均滤波已经在图像数据、雷达数据、 声音数据的随机噪声压制中得到了广泛的应用， 但该 方法目前还未见在叠前地震数据去噪中得到应用。 ChaoXing 88 煤田地质与勘探 第 42 卷 原始的非局部平均算法计算每一个点滤波系数 时要对整个数据体进行权值运算，因此这种方法的 计算量非常大，难以适用于叠前地震资料处理。为 使得非局部平均滤波应用于叠前处理，需要减少权 的计算次数。在图像处理领域，研究人员提出了许 多提高计算效率的方法。 Mahmoudi 等[2]利用图像子 块的均值与平均梯度方向对邻域进行预选择，生成 相似像素集合，这种方法可以显著提高运算速度， 但是对于弱梯度与低信噪比图像去噪效果不佳； Coupe 等[3]采用图像子块的均值与方差选择与当前 图像子块相似的邻域，其对于低信噪比资料有较好 的适应性，但存在改变纹理走向的风险；Ramana- than[4]利用概率统计的方法排除不相似子块， 该方法 对于低信噪比资料适应性较差； Tasdizen 等[5]利用主 成份分析将特征向量空间用于相似系数求取，这种 矩阵向量化的方法可能破坏图像的结构信息；Brox 等[6]利用聚类树生成相似集，该方法存在存储代价 过高的问题。 本文通过分析叠前地震数据的结构冗余特性， 对原始非局部平均算法进行了改进，提出了可以适 用于叠前地震数据处理的方法。这主要包括叠前地 震数据体分块、相似邻域预处理、相似邻域优选、 自适应滤波参数选择 4 部分内容。 1 非局部平均滤波基本原理 设原始地震记录为f，由有效波 s 与干扰波 n 组成，即fsn 。 NMLF 处理后的记录为 ˆ f，i、 j 为地震记录中的采样点，点i的滤波后结果为若干 与其局部结构相似的样点j加权求和得到。每一点 j均具有特定的权系数。 为定量计算i与j点的结构 相似度，需要定义一个包含结构信息的邻域作为比 较单位。邻域的维度由数据的维度决定，其形状和 大小是可变的。i点与j点的局部结构相似度是利用 邻域间高斯加权的欧几里德距离 , D i j来表示。 设以i点为中心的邻域为 N i，滤波后样点的 振幅为 ˆ f i，以j点为中心的邻域为 N j，滤波前 j点的振幅为 f j， , w i j为权系数， , D i j为高 斯加权欧式距离，A为搜索窗口，h为滤波参数， a G为高斯核函数，a 为高斯核函数的标准差， 0 x、 0 y表示高斯核的中心点坐标。将 NMLF 用公式表 示如下 ˆ , j A f iw i j f j 1 , , exp D i j w i j h - 2 2 2 , aij Di jGf Nf N- 3 22 00 , exp 2 a xxyy Gx y a -- - 4 其中，滤波参数h控制指数函数的衰减速度，起着 控制平滑程度的作用。 高斯核函数确定i的邻域与j 的邻域的欧几里德距离，给予远距离的点较小的权 值，近距离点的邻域赋以较大的权值，使得与中心 点邻域最接近的结构保存表示集合内对应元素 相乘。 , w i j满足0 , 1w i j≤≤并且 , 1 j w i j 。 其计算方式如图 1 所示，五星所在点为计算点采样 点，方框内的部分为邻域，3 种不同的线段代表 3 种不同的局部结构。对于任意一个计算点p，其局 部结构可用以其为中心的邻域表示，滤波时与该点 邻域结构越相似，则分配的权值越大，局部结构相 同的计算点权重最大。图中p点与1q点的邻域结构 相似度最高，因此其权值 , 1w p q远大于2q的权值 , 2w p q与3q的权值 , 3w p q。这样，参与滤波点 的权重与其位置无关，仅仅与计算点邻域的局部结 构相似性有关。因此，这种滤波方法称为非局部平 均滤波算法。NLMF 算法正是利用有效信号局部结 构重复出现的概率大于干扰信号来增强有效波压制 干扰波。 图 1 NMLF 加权计算示意图 Fig.1 The weighted of NMLF 原始的 NMLF 算法中搜索窗口A为整个数据 体。设数据体、搜索窗口、邻域的大小分别为 x、y 和 z，那么 NMLF 的计算复杂度为O xyz，因此选 择合适的搜索窗口与邻域可以显著提高计算效率。 2 改进的 NMLF 方法 2.1 搜索窗口的选择 非局部平均滤波计算中计算量最大的是权系数 的计算，因此仅选择对滤波效果影响较大的邻域进 行权系数计算是一种提高算法效率的方法。但这种 方法是否有效取决于数据的预选择特征与采用的预 选择策略。相比图像资料，叠前地震数据的噪音背 ChaoXing 第 5 期 胡新海等 叠前非局部平均滤波压制随机噪音 89 景要大得多，但其规律性较强，结构复杂度较低。 基于这些特点， 在 CMP 道集中一段偏移距内开一个 时窗作为运算邻域，然后利用块匹配规则在搜索窗 口中找到该 CMP 道集与相邻 CMP 道集中与该邻域 近似度满足某个阈值的邻域，再将这些邻域堆叠成 一个三维矩阵；用可分的三维变换系数矩阵表示该 三维矩阵，利用收缩变换系数进行去噪处理，然后 再进行三维逆变换，最终将所有数据块与系数进行 聚合得到去噪后的结果。 对于叠前 CMP 道集， 其搜索窗口为一系列时空 窗的叠合图 2。首先进行数据预处理，将叠前地震 数据按照 CMP 线号、CMP、偏移距进行分选，对 CMP 道集做动校正；然后根据速度谱上控制标志层 的谱点进行时间窗划分，其标准是保证主要目的层 结构特征相近的部分都处于同一个搜索窗内。图 2 中速度谱上划分的时间窗包含 CMP1、 CMP2、 CMP3 三个 CMP 点， 分别对应 3 个道集时间方向的窗口范 围；在偏移距方向可以选择全偏移距也可以对偏移 距进行等分。依照同样的办法，将 n 个相邻面元的 CMP 道集进行相似的划分，然后将得到的搜索窗口 进行叠合，形成三维搜索域。其中 n 的大小取决于 水平方向的构造平缓程度。需要注意的是，当构造变 化较大时，相邻 CMP 道集同一位置相似性较差，n 的 值过大会使去噪效果变差，并可能出现结构噪声。 图 2 搜索窗口选择示意图 Fig.2 The for selection of scan windows 2.2 相似集的选取 相似邻域的选择依赖于对数据局部特征的提 取。 经过静校正与动校正之后的叠前 CMP 道集有效 信息近似于水平层状，利用奇异值分解算法难以表 示因为 AVO 效应或构造因素引起的振幅变化。因 此，选择梯度域奇异值分解法进行有效信息表示。 对于搜索窗口内的数据，首先将三维搜索域中的各 个搜索窗口 W 分别计算数据的梯度，然后将其排列 成2N的矩阵 D 进行奇异值分解。 [1 ,2 ,3 ,, ] TTT T TT DWWW W NUV ⋅⋅⋅⋅  5 其中 [,]T W iW i W i xy  为 i 点的梯度；U是 NN的正交矩阵；是2N的奇异值矩阵；V为 大小22的正交矩阵，其第二列向量对应最小的奇 异值。奇异值矩阵的特征向量 1 λ、 2 λ反映特征向量 方向的能量变化。在有效波位置振幅变化较小，特 征值近似为零；在“断点”边缘位置振幅变化大，特 征值 12 λλ0。 由于奇异值具有较强的抗噪能力， 且噪声没有任何方向性，因此 1 λ、 2 λ能够表征数据 的 局 部 特 征 。 采 用 文 献 [7]中 提 出 的 方 法 ， 将 12 λ λ e iii作为局部方向能量的度量， 并进行 归一化处理 min max min ee e ee - - 6 e表示图像的局部结构，e较大的区域包含丰富的 细节信息，e较小的区域相对比较平坦。将e作为 预选择的依据，计算各图像子块的e生成结构特征 图。 对于每一点， 首先比较当前邻域的e与其他点e 的差值，生成相似集[7-8]。相似集的大小与邻域的形 状大小是影响滤波的效果的重要因素，需要进行试 验来确定。 2.3 自适应滤波参数选择 不同 CMP 道集的信噪比不同，即使同一 CMP 道集内不同位置的信噪比也有高有低。当信噪比水 平超过一个数值以后， 不可能找到一个全局的滤波 参数使其对图像的各个部分都具有很好的去噪效 果[9-10]。滤波参数h直接决定 NLMF 的滤波性能， ChaoXing 90 煤田地质与勘探 第 42 卷 对h进行自适应处理使其能够随信噪比变化而变化 以达到好的去噪效果。根据对地震数据的冗余度分 析可以得出信噪比高的成分结构冗余度高，信噪 比低的成分结构冗余度低。因此，可以根据数据冗 余度代替信噪比来控制滤波参数的大小。由权值表 达式 , , exp D i j w i j h -可以看出 邻域之间的欧 式距离与滤波参数的比值决定了权值的大小。冗余 度不同的点，其邻域之间的欧式距离具有不同的分 布，为了得到更好的滤波效果，对于不同的距离分 布应该选用不同的滤波参数。根据冗余度将数据分 为冗余度大的邻域和冗余度小邻域。以数据格点邻 域之间的欧式距离小于某个阈值的个数 num 占搜索 窗点数win的比例p作为分类的准则，即 pnum win 7 则对于不同冗余度的邻域，有 max λhpd a ⋅⋅ 8 式中 h为滤波参数；λ为基础加权系数；a为衰减 指数； max d为相似集中最大的欧几里德距离。 资料中冗余度高的邻域p值较大，对应较大的滤 波参数。a一般为正值，a越大则不同冗余度的邻域 间滤波参数大小差别越大。欧几里德距离越大代表冗 余度越高，不同位置的最大欧几里德距离不同，可以 利用控制邻域的 max d对应的权值来控制整体权值的分 配。式8综合了滤波效果的整体调整和局部调整，利 用λ和a来控制整体效果，同时利用不同位置p与 max d对于不同位置不同结构的邻域自适应地选择滤波 参数，而不是全局利用一个固定参数，这样可以避免 有的位置过度平滑而有些位置滤波不充分的情况。 4 应用效果 图 3 为模型数据滤波效果图，由此表明相比 常用的 FX 预测反褶积RNA， 非局部平均滤波对于 随机噪声具有更加有效的压制作用。 图 3 模型数据非局部滤波效果图 Fig.3 Effect of nonlocal filtering of model data ChaoXing 第 5 期 胡新海等 叠前非局部平均滤波压制随机噪音 91 图 4 为某区块的实际地震资料。图 4a 为动校正 后未做非局部平均滤波的道集， 图 4b 为非局部平均滤 波后的道集，可见非局部平均滤波可以较好的压制随 机噪音，显著提高地震数据信噪比，增强有效波的同 相性；图 4c 为保幅性说明图，实线为原始道集 a 时窗 范围 2 7002 800 ms 的 AVO 曲线， 虚线为滤波后道集 b 时窗范围 2 7002 800 ms 的 AVO 曲线。滤波前后 a 道集的 AVO 曲线有所改变，而相对振幅关系得到了 保持，由此可见非局部均值滤波具有保幅性。图 4d 原始道集的叠加剖面， 图 4e 为非局部平均滤波后道集 的叠加剖面，表明经过叠前非局部平均滤波处理后， 剖面整体信噪比明显提高，层位连续性增强。 图 4 实际资料非局部滤波效果图 Fig. 4 The result of NMLF on actual seismic data 5 结 语 提出的用于叠前随机噪音压制的快速非局部平 均算法，采用梯度域奇异值分解提取数据的局部结 构信息，构造局部结构特征描述因子，在此基础上， 提出了相似集优选方法，求取体现结构相似度的权 系数，并结合地震数据的特点，利用速度场自动划 分搜索窗口；通过算法分析，提出了自适应滤波参 数的选择方法代替原有的单一参数，能够对数据不 同信噪比的局部区域更好地进行滤波；应用模型地 震数据与实际地震资料进行测试，验证了方法的保 幅性与有效性。 该方法适用于低信噪比地区的三维叠前地震数 据区处理，可以提高剖面整体信噪比，增强层位连 续性。叠前地震资料数据量大，而且该方法属于全 局搜索算法，虽然进行了优化，但运算与存储量依 然较大，因此要将算法进行并行化改造作为努力的 方向。 参考文献 [1] 崔树果，朱凌燕，王建花. f-x 域 Cadzow 技术分块压制随机噪 声及其应用[J]. 石油物探，2012，51144–49. [2] MAHMOUDI M，SAPIRO G. Fast image and video denoising via nonlocal means of similar neighborhoods[J]. IEEE Signal Processing Letters，2005，1212839–842. [3] COUPE P，YGER P，PRIMA S，et al. An optimized block- wise nonlocal means denoising filter for 3-D magnetic images[J]. IEEE Transactions on Medical Imaging，2008， 274425–441. [4] VIGNESH R， OH B T， KUO C C J. Fest non-local meansNLM computation with probabilistic early termination[J]. IEEE Signal Processing Letters，2010，173277–280. [5] TASDIZEN T. Principal neighborhood dictionaries for nonlocal means image denoising[J]. IEEE Transactions on Image Processing，2009，18122649–2660. [6] BROX T，KLEINSCHMIDT O，CREMERS D. Efficient nonl- ocal means for denoising of textural patterns[J]. IEEE Transa- ctions on Image Processing，2008，1771083–1092. [7] 许光宇， 谭结庆， 钟金琴. 有效保持细节特征的快速非局部滤 波方法[J]. 计算机工程与应用，2012，4823196–202. [8] 余松煜. 数字图像处理[M]. 北京电子工业出版社，1989 18–57. [9] 许录平. 数字图像处理[M]. 北京科学出版社，2007. [10] 刘涛. 小波域中的非局部平均去噪算法研究[D]. 西安 西安电 子科技大学，2010. ChaoXing