中值滤波是一种非线性滤波方法^[1]，并广泛的应用于地震资料噪音压制^[2]和上下行波波场分离^[3-4]等。中值滤波压制脉冲噪声具有明显优势，但同时造成有效信号能量的损失^[5]。针对中值滤波方法存在的问题，刘洋等^[6]提出一维时变中值滤波，优化了脉冲噪声的压制效果。王伟等^[7]为了适应地震数据的时空变化，对中值滤波窗口长度进行自适应调整，实现随机噪音压制；周丽等^[8]利用多方向中值滤波迭代计算实现同步激发震源波场分离；于富文等^[9]利用多方位矢量中值滤波对可控震源同步激发采集干扰进行压制。为了提高压噪能力和减少对有效信号的损伤，上述研究在常规中值滤波的滤波窗口长度和滤波方向等方面进行优化。

开关非局部中值滤波^[10](Switching Non-local Median Filter)考虑到地震数据结构的相似性，借鉴非局部均值的处理思路，避免了均值滤波对有效信号的平滑效应。开关^[11](Switch)是指通过噪音点识别，有针对性进行噪音处理，从而避免对非噪音点的过度损伤；非局部中值滤波^[12](Non-local Median Filter)在计算目标点的数据块与其周围临近数据块的相似度的基础上，按照相似度的大小对临近数据块的中心点进行排序，选择一部分临近数据块的中心点作为候选点，使用候选点的中值代替待处理点的值，完成噪音的处理。本文针对海洋浅地层剖面模拟数据与实际资料，着重分析讨论各参数对方法的影响，以达到优选参数组合压制随机噪声的目的。

1 方法原理

开关非局部中值滤波^[10](SNLM)利用二维数据结构的相似性进行随机噪音的去除，包括噪声检测和噪声去除两部分，噪声检测包括异常点检测和梯度值检测。首先通过异常点检测对数据部分的噪音点进行筛选，而对于异常点检测未响应的噪音点，利用梯度值检测进行目标数据点的孤立性判断。最终将满足检测条件的目标数据点判断为噪音。X、Y、Z分别为原始数据、非局部中值滤波器输出的数据、噪音去除后的结果(见图 1)。

1.1 非局部中值滤波

SNLM包括非局部中值滤波(NLM)和噪音点检测(Switch)两部分，由于噪音点检测需要用到NLM的计算结果，因此首先介绍NLM的计算步骤^[12]：

(1) 假设目标点为p(m, n)(m、n是p所在数据中的位置)，以p为中心，划定一个宽度为r的正方形区域，记为S，即目标数据块；再以p为中心，划定一个宽度为ρ的正方形区域，记为I，即搜索数据块(相似度计算区域)，其中宽度r与宽度ρ均为奇数，且ρ>r。

(2) 分别计算目标数据块S和对比数据块的相似度，计算过程按照对应点差值的平方和作为相似度的评价指标，计算公式如式(1)。其中对比数据块与S维度相同，且其中心点是搜索数据块I中的所有的点。假设所有对比数据块中心点的值依次为a_{1, c}, a_{2, c}, …, a_{ρ²－1, c}。

$ e_j=\sum\limits_{k=1}^{r^2-1}\left(x_k-a_{j, k}\right)^2, j=1, 2, \cdots, \rho^2-1 。$

(1)

式中：x_k为目标数据块S内除p(m, n)之外的其它数据点的值；a_{j, k}为第j个对比数据块内除其中心点的其它数据点的值；e_j为计算的相似度值。

(3) 按照第(2)步的计算，得到ρ²－1个相似度值。按照从小到大排序后的相似度值对应的中心点值的顺序变为a_{1, c}^*, a_{2, c}^*…, a_{ρ²－1, c}^*，并选择数组前N个数据点作为中值滤波的候选点，即a_{1, c}^*, a_{2, c}^*…, a_{N, c}^*，将候选点的中间三个点的均值作为NLM的输出y，即：

$ y=\operatorname{mean}\left(a_{\frac{N-1}{2}, \mathrm{c}}^*, a_{\frac{N+1}{2}, \mathrm{c}}^*, a_{\frac{N+3}{2}, \mathrm{c}}^*\right) 。$

(2)

1.2 噪音点检测

SNLM的第二个环节就是噪音点检测。噪音点检测^[11]包括两部分：(1)计算目标点和NLM输出的相对误差值；(2)目标点的孤立性检测。首先假设目标点的值为p_{m, n}，NLM输出的数据点值为y_i，噪音点检测公式如式(3)。

$ f_{\mathrm{p}}=\left\{\begin{array}{l} 1, \left|p_{m, n}-y_i\right| /\left|y_i\right| \geqslant e r r \\ 0, \left|p_{m, n}-y_i\right| /\left|y_i\right| <e r r \end{array}\right. \text {。} $

(3)

式中：err是门限阈值；y_i是NML输出结果；f_p是检测结果。因为地震有效反射信号与随机噪声的振幅值差别较大，所以将err设置为比例系数。当目标数据比值超过设定的阈值时，判定目标点是噪音点；而目标数据比值小于阈值时，需要进一步验证目标数据的孤立性，若孤立性成立，亦判定为噪音点，若孤立性不成立，该目标数据为有效信号。

孤立性检测首先在以目标点为中心点的3×3的矩形区域内，计算目标点不同方向的梯度值，再对其各个方向的梯度变化进行检测(式4至式8)。

$ g_{n, \mathrm{~h}}=\left\{\begin{array}{l} 1, \left(p_{m-1, n}-p_{m, n}\right)\left(p_{m, n}-p_{m+1, n}\right)<0 \\ 0, \text { 其他 } \end{array}, \right. $

(4)

$ g_{m, \mathrm{v}}=\left\{\begin{array}{l} 1, \left(p_{m, n-1}-p_{m, n}\right)\left(p_{m, n}-p_{m, n+1}\right)<0 \\ 0, \text { 其他 } \end{array}\right. \text {, } $

(5)

$ g_{d 1}=\left\{\begin{array}{l} 1, \left(p_{m-1, n-1}-p_{m, n}\right)\left(p_{m, n}-p_{m+1, n+1}\right)<0 \\ 0, \text { 其他 } \end{array}\right. \text {, } $

(6)

$ g_{d 2}=\left\{\begin{array}{l} 1, \left(p_{m-1, n+1}-p_{m, n}\right)\left(p_{m, n}-p_{m+1, n-1}\right)<0 \\ 0, \text { 其他 } \end{array}, \right. $

(7)

$ g_{\mathrm{p}}=g_{n, \mathrm{~h}}\left|g_{m, \mathrm{v}}\right| g_{d 1} \mid g_{d 2} 。$

(8)

式中：(m, n)是目标点所在的坐标；p代表各个点上的数值；g_{n, v}是垂直方向的检测结果；g_{m, h}是水平方向的检测结果；g_d1和g_d2是倾斜方向的检测结果；g_p是目标点的孤立性最终检测结果；如果g_p=1则判定目标点是孤立点。

综上分析，异常点检测与梯度值检测的判定值满足至少一项为“真”时，目标点数据将被认定为噪音点，反之，目标数据被判定为有效信号。即该噪音点检测方法可以表示为：

(1) 目标点与NLM输出的差值大于等于门限阈值。

(2) 目标点被判断为孤立点。

$ Z_{m, n}=\left\{\begin{array}{l} y_i, f_{\mathrm{p}}=1 \text { 或 } g_{\mathrm{p}}=1 \\ p_{m, n}, \text { 其他 } \end{array}。\right. $

(9)

式中Z_{m, n}是滤波结果。

2 关键参数分析

滤波参数的选取决定剖面随机噪声的压制效果。为达到最优的滤波效果，建立理论模型，选择主频为300 Hz的零相位雷克子波，采样间隔0.1 ms，1个曲线型同相轴、4个直线型同相轴。其中倾斜同相轴和曲线型同相轴右侧的反射系数为负值，其余同相轴反射系数为正值，下层的水平同相轴反射系数从左至右渐变(见图 2a)。在剖面内加入随机噪声(见图 2b)。分析该方法的关键参数目标数据块宽度r、搜索数据块宽度ρ、相似度最优的数据点数N、噪音检测的阈值参数err对噪声压制效果的影响。采用频域SVD方法^[13]对去噪后剖面进行信噪比估算，分析不同参数对噪声压制效果；计算有效信号与压噪后剖面的互相关系数corr1，反映噪音压制前后的波形相关性；计算有效信号与去除的噪声剖面的互相关系数corr2，反映滤波前后有效信号的损伤。通过分析估算信噪比和两种互相关系数的变化趋势，讨论四个关键参数的适用范围。

2.1 目标数据块宽度r

因为必须满足ρ>r的条件，目标数据块宽度必须先确定。取参数r为3~17范围内的奇数，其余三个参数分别设定为ρ=21，N=19，err=0.2。图 3为选取滤波参数r为3、9和17的滤波结果。从图 3中红框可以看出，当r值取9时，噪声压制效果优于其他，尤其以最上方的横向同相轴能量变化明显。

分析去噪后剖面的估算信噪比、有效信号与去噪后剖面的相关系数1、有效信号与去除的噪声剖面相关系数2等随着r值增大的变化趋势(见图 4)。r值越大，估值信噪比与相关系数1均增大，而相关系数2减小。当r值增加超过13时，相关系数2不降反增，即去除的噪声剖面与有效信号的相关性增加，这说明有效信号的能量损失。上述分析结果表明，r值过小，目标数据点的特征无法充分表现，噪音压制不充分；而当r值接近搜索数据块区域宽度ρ时，目标数据点的局域性降低，数据聚焦不充分，造成相似度计算产生偏差，进而滤波结果不准确且损失有效信号能量(见图 3c中蓝框)。针对模型数据选择目标点宽度r值为11。

2.2 搜索数据块宽度ρ

搜索数据块宽度ρ的实质就是窗口宽度，优选合适的搜索数据块宽度不仅影响运算的效率，更影响噪音压制效果。设定其余三个参数依次为r=11，N=19，err=0.2。分别计算ρ为15~41范围内部分奇数的滤波结果，图 5选取ρ为15、19和29滤波结果。当ρ值偏小时，纯噪音区域随机噪音残留严重，同相轴能量恢复不明显。ρ值逐渐增大时，滤波效果也随之变好，但增加到一定程度后，滤波效果改善不大。分析去噪后剖面的估算信噪比、有效信号与去噪后剖面的相关系数1、有效信号与去除的噪声剖面相关系数2等随着ρ值增大的变化趋势(见图 6)。ρ值越大，估值信噪比与相关系数1均增大，而相关系数2减小。但在ρ值大于21后，估算信噪比、相关系数1和相关系数2的变化趋于平缓。这说明搜索数据块宽度增大至一定程度后，数据相似度将不再改善。并且随着ρ值的增大，每个时间点计算相似度的次数增大，计算时间增多。因此在选择参数ρ时，需要兼顾噪音压制效果和计算效率。针对模型数据选择参数ρ值为21。

2.3 相似度最优的数据点数N

当选择合适的相似度计算区域大小ρ和目标数据点所在数据块的大小r时，作为判断参数相似度最优的数据点数N影响中值滤波中的平滑效应。设定其余3个参数分别为r=11, ρ=21，err=0.2，分别计算N在5~31范围内奇数的滤波结果，图 7是N为5、19和31的滤波结果。在其余参数都一致的情况下，相似度最优的数据点数较小时，噪音压制效果不佳，其原因是中值滤波候选点数有限，即使计算出足够多的和目标数据点所在数据块相似度大的数据块，仍然削弱了中值滤波的效果。而如果N选择偏大，必然会将中值滤波的平滑效应放大，存在的模糊效应增加，损伤有效信号。随着N的增大，估算信噪比、相关系数1与相关系数2均呈现增大趋势。N值增大到19之后，估算信噪比与相关系数1变化趋势放缓，而相关系数2仍然较快增大(见图 8)。为避免N值增大带来的模糊效应，针对模型数据选择参数N值为19。

2.4 噪音检测的阈值参数

在噪声检测中使用的阈值参数，其大小制约检测出的噪音边界。设定其余三个参数分别为r=11，ρ=21，N=19，分别计算err在0.1~0.5范围内的滤波结果。阈值参数err是利用振幅比值区分随机信号与有效信号。原始剖面的信噪比影响阈值参数的选择，原始剖面中随机噪声振幅比例大，则选择较大的阈值参数，反之选择较小的阈值参数。与其他三个参数来比较，随着阈值参数的增大，估值信噪比、相关系数1与相关系数2均呈现下降趋势，并且变化范围有限。针对模型数据，阈值参数取0.2时，信噪比略高(见图 9)。实际数据处理过程中，应当根据实际资料的信噪比，进行阈值参数的优选。信噪比小对应选择较大阈值参数。反之，选择较小阈值参数。

综合上述分析可知，目标数据块宽度r与搜索数据块宽度ρ相互制约，r值决定数据块特征是否明显，ρ值决定搜索数据是否足够。当r值远小于ρ值，目标数据特征不足，无法对噪声信号进行识别，滤波效果不好；当r值接近于ρ值，搜索数据量不足，有效信号损失明显(见图 3c与图 5c中的蓝框)。从参数分析结果中得出，r值与ρ值的比值略高于50%，滤波效果最优。在数据块边界搜索过程中，会有能量的泄漏。边界数据的计算应将已处理样点的滤波结果作为搜索区域中的样点，压制随机噪音的同时保护有效信号。

3 应用结果

针对理论模型，分析目标数据块宽度r、搜索数据块宽度ρ、相似度最优的数据点数N、以及噪音检测的阈值参数err等关键参数对随机噪声压制的效果，数据信噪比的改善，以及有效信号的损伤等方面的影响，选择滤波参数分别为r=11，ρ=21，N=19，err=0.2。从图 10a可以看出，SNLM能够有效滤除随机噪音，同相轴的连续性和能量改善，剖面信噪比得到很好的提升。

选择原始未经压制随机噪声的浅地层剖面实际数据，进行方法有效性验证。该数据是在渤海近岸的河流入海口附近采集的，环境噪音比较强。根据参数优选原则，保持目标数据块宽度与搜索数据块宽度值的比例，相似度最优的数据点数与阈值参数交互选取，以逐步多次的方式进行实际数据的滤波参数选取，确定r=13, ρ=27，N=21, err=0.2的滤波参数。图 11是原始剖面和去噪后的剖面。从图 11中可以看出，原始数据受噪音的影响比较严重，局部同相轴模糊不清。

利用本文滤波方法处理后，信噪比提高，同相轴连续性变好，随机噪声压制效果较好。对比处理前后的数据剖面，小能量的反射同相轴没有在滤波后被模糊化，同相轴断点位置也更加清晰。图 11c的噪音剖面中随机噪声特征明显，无明显的有效信号损失。

4 结论

本文应用开关非局部中值滤波对浅地层剖面资料进行随机噪声压制，获得了较好的效果。建立模型数据，着重分析了目标数据块宽度r、搜索数据块宽度ρ、相似度最优的数据点数N、以及噪音检测的阈值参数err等关键参数，优选出了针对模型数据的滤波参数，并总结出关键滤波参数的选择原则：

(1) 目标数据块宽度r与搜索数据块宽度ρ相互制约。r值太小会导致目标数据特征拾取不足；r值过大接近于ρ值会导致搜索数据量不足，损失有效信号。r值与ρ值的比例略高于50%最佳。

(2) 相似度最优的数据点数N、噪音检测的阈值参数err与数据信噪比相关，应根据信噪比大小进行上述参数的选择。

(3) 搜索数据块边界计算问题能够导致信号能量的泄露，参照临近样点滤波结果进行搜索，不仅能够更好的保护有效信号，而且压制随机噪音效果更好。