2019, Vol. 62
2. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
2. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum-Beijing, Beijing 102249, China
地震资料在采集过程中由于人为、环境以及仪器的影响会引入随机噪声，它会降低地震资料的信噪比，影响后续的处理及解释工作，因此希望通过一定的处理尽可能地消除掉随机噪声，以提高信号估计的精度.实际生产中，常用的压制噪声的方法有F-X反褶积(Gulunay, 1986)、多项式拟合(Yu et al., 1989)、F-K滤波(张军华等, 2005)、奇异值分解(刘伟等, 2010)和小波变换(王勇等, 1998)等，然而这些方法在去噪处理过程中都存在不同程度的局限性，尤其在复杂构造及低信噪比情况下压制随机噪声的同时会损失一部分有效信号或产生假象，降低去噪剖面的保真度.随机信号与有效信号的最大差别就是是否具有连续性，随机信号在时空中杂乱分布，而有效信号则以一定的规律存在即我们所说的同相轴，所以我们可以将有效信号同相轴看作是一种曲线，从而利用曲线检测算法分离有效信号与随机噪声.考虑到Hessian矩阵(Frangi et al., 2000)在检测曲线结构上表现出的良好性能(Koller et al., 1995)，本文提出了一种基于Hessian矩阵的多尺度信号增强算法.对于图像某点的Hessian矩阵，其最大特征值和其对应的特征向量对应其邻域二维曲线最大曲率的强度和方向，最小特征值对应的特征向量对应与其垂直的方向，简单来说，图像某点的Hessian矩阵特征值大小和符号决定了该点邻域内的几何结构.由此原理设计的一些滤波器在医学领域的MRA和CTA血管造影增强中有很重要的应用(Lorenz et al., 1997; Xu et al., 2013; Shashank et al., 2013).我们可将地震数据中的同相轴看作为一种特殊的“血管”，从而利用Hessian矩阵的特性增强我们的有效信号，达到压制随机噪声的目的，并且由于该方法不受同相轴倾角变化的限制，所以在去除随机噪声的同时能较好的保护复杂构造中的有效信号.从本文中设计的模型及实际资料处理结果来看，本方法取得了较好的去噪效果，与传统去噪结果的对比也说明了新方法在噪声压制及有效成分保持方面的优势.
1 Hessian矩阵的应用原理Hessian方法是一种利用高阶微分提取图像特征的方法(游嘉和陈波, 2011).分析图像L局部特征的常用方法是在点x0的邻域中考虑它的泰勒级数展开，即
|
(1) |
式中，Δx0为点x0在其某个邻域内的变化量；
|
(2) |
H(x0)中的元素
|
(3) |
以及
|
(4) |
通过分析式(3)和(4)我们可以得到如下几何解释，由Hessian矩阵的特征值分解得到两个正交的方向向量，且在一定尺度因子下保持不变.假设λ1和λ2是两个特征值，ν1和ν2是对应的两个特征向量，
|
表 1 Hessian矩阵特征值和探测形状之间的对应关系 Table 1 Correspondence between eigenvalues and detected shapes of Hessian matrix |
在这里我们采用Frangi的增强滤波函数：
|
(5) |
其中P(σ, x0)是点x0的响应度，Rb=λ1/λ2，β和c分别是控制线性特征及整体平滑的参数.
2 基于Hessian矩阵的多尺度噪声压制地震数据中含有不同频率的信号，也就是同相轴的粗细程度不同，因此我们需要在多个尺度上进行线性增强滤波以达到去噪的需求.将Hessian矩阵的差分运算与高斯函数结合，通过改变高斯函数的标准偏移量就可以获得不同尺度下的线性增强滤波.根据线性尺度空间理论(Koenderink, 1984; Florack et al., 1992)，尺度空间函数由输入图像与高斯滤波器的二阶导的卷积得到：
|
(6) |
高斯函数G(x, y, σ)表达式为：
|
(7) |
获得不同的尺度空间函数后分别计算在各个尺度上的滤波器响应结果，最终的响应函数可由下式确定：
|
(8) |
为了进一步说明多尺度理论在基于Hessian矩阵去噪中的必要性，设计了如图 1所示的数值模型.图 1a为设计的数值模型及其频谱，包含了三条水平同相轴，从上至下主频分别为70、50和30 Hz.1b—1e是尺度因子分别等于2、6、10、14时的高斯滤波器及线性响应输出结果，图 2是各个尺度下线性响应输出结果对应的频谱.可以看出，当尺度因子较小时，线性滤波器响应输出的主频较高；尺度因子增大时，线性滤波器响应输出的主频随之降低.由此我们可以得出如下结论：基于Hessian矩阵的多尺度滤波器在小尺度时输出高频信号，大尺度时输出低频信号，因此，我们对地震资料进行去噪处理时，要想兼顾有效信号高频和低频的能量就必须在多个尺度上进行处理，获得不同的尺度空间中的滤波器响应，最终综合得到去噪结果.
|
图 1 Hessian矩阵多尺度响应结果测试. (a)数值模型及其频谱; (b)—(e)分别为尺度因子σ为2、6、10、14时的高斯滤波器及线性响应输出. Fig. 1 Test of Hessian matrix multi-scale corresponding result. (a) Numerical model and its spectrum; (b)—(e) Gaussian filters and linear response with σ equal to 2, 6, 10 and 14 respectively. |
|
图 2 各尺度下的线性响应输出对应的频谱 Fig. 2 Linear response output spectrum corresponding to each scale |
因此，基于Hessian矩阵的多尺度随机噪声压制方法实现步骤可概括如下：
(1) 输入原始数据L，设定尺度因子σ范围
(2) 初始化尺度因子和响应函数P(σ, x0)max；
(3) 计算元素Lij与高斯函数二阶微分的卷积，生成Hessian矩阵，并计算其特征值和特征向量；
(4) 计算响应函数P(σ, x0)；
(5) 迭代尺度因子σ=σ+Δσ，转至步骤(3)；
(6) 迭代结束，输出响应函数最大值P(σ, x0)max；
计算振幅补偿因子A，得到最终的滤波结果
其中，A是一常数，可通过处理前后剖面的平均值计算得到，如式(9)所示：
|
(9) |
为了验证本方法的去噪性能，设计了如下数值模型，其中图 3a为原始未加噪声数据，共390道，500个采样点，采样间隔1 ms.图 3b为加入带限高斯随机噪声数据(噪声频宽与有效信号频宽一致)，信噪比为-0.0093 dB.分别采用了本文提出算法和传统算法F-X反褶积进行了去噪处理，结果如图 4所示.
|
图 3 (a) 无噪数据; (b)含噪声数据 Fig. 3 (a) clean data; (b) noisy data |
图 4a、b和c分别是基于Hessian矩阵的单尺度、多尺度和F-X反褶积去噪结果.其中，基于Hessian矩阵的单尺度和多尺度去噪结果的信噪比分别为3.7356 dB和8.1209 dB，F-X反褶积去噪结果的信噪比为6.3408 dB.可以看出，两种方法对噪声都有一定的压制效果，基于Hessian矩阵的单尺度去噪效果明显不如多尺度去噪效果，单尺度去噪只保留了某一频率附近的有效信号，差异剖面中可以明显的看到同相轴的存在.而多尺度去噪则解决了这一问题，绝大部分的有效信号被有效的保留了下来，差异剖面中几乎看不到有效信号的存在.F-X反褶积算法对复杂构造中的有效信号能量不能起到很好的保护作用(如图 4f中箭头所示)，并且F-X反褶积去噪后产生了一些虚假的能量较弱的同相轴，降低了去噪剖面的保真度，这对后期进一步的处理解释也带来了影响.总体来看，基于Hessian矩阵的单尺度的去噪算法不能保护所有频率的有效信号，多尺度的去噪算法则较好的解决了这一问题，并且与F-X反褶积算法相比，基于Hessian矩阵的多尺度去噪算法的去噪结果不仅信噪比更高，对随机信号压制更为彻底，而且保护了复杂构造中的有效信号，剖面的保真度也较高.
|
图 4 本文所提算法的去噪结果和F-X反褶积去噪结果对比. (a)基于Hessian矩阵的单尺度去噪结果；(b)基于Hessian矩阵的多尺度去噪结果；(c) F-X反褶积去噪结果；(d)基于Hessian矩阵的单尺度去噪的噪声数据；(e)基于Hessian矩阵的多尺度去噪的噪声数据；(f) F-X反褶积去噪的噪声数据 Fig. 4 Comparison of denoising results by the proposed method and F-X deconvolution. (a) Single-scale denoising results based on Hessian matrix; (b) Multi-scale denoising results based on Hessian matrix; (c) Denoising results of F-X deconvolution; (d) Noisy data of single scale denoising based on Hessian matrix; (e) Noisy data of multi-scale denoising based on Hessian matrix; (f) Noisy data using F-X deconvolution. |
为了进一步检验基于Hessian矩阵的多尺度去噪算法的性能，选取了某工区实际地震资料进行处理，如图 5所示，数据共396道，1001个采样点，采样间隔2 ms.由于随机噪声的存在，使得同相轴连续性较差，信噪比整体也较低，进一步影响了资料的分辨率，并且该资料的构造比较复杂，尤其是中深层大斜率同相轴的存在给资料处理带来了较大的难度.同样，我们对此剖面利用本文所提算法和F-X反褶积算法进行了去噪处理，结果如图 6所示.
|
图 5 原始地震数据 Fig. 5 Original seismic data |
|
图 6 叠后地震资料去噪结果对比. (a)本文算法去噪结果；(b) F-X反褶积去噪结果；(c)本文算法去噪后的噪声；(d) F-X反褶积去噪后的噪声 Fig. 6 Comparison of denoising results of post stack seismic data. (a) Denoising results of the proposed algorithm; (b) denoising results of F-X deconvolution; (c) Noisy section using the proposed method; (d) Noisy section using F-X deconvolution |
图 6a和b分别是基于Hessian矩阵的多尺度去噪和F-X反褶积去噪结果，图 6c和d是两种去噪算法得到的噪声剖面.从实际去噪效果来看，F-X反褶积算法虽然具有一定的去噪效果，但是对有效信号有不同程度的损伤，从噪声剖面中可以明显看到有效信号的存在(如6d中红色箭头所示).而本文方法去噪结果信噪比更高，而且有效信息也更丰富，同相轴相对更加连续，噪声剖面中基本没有有效同相轴的痕迹，说明该方法对有效信号的损伤较小.图 7为两种方法去噪后的频谱对比，可以看到，本文所提算法去噪结果的频谱与F-X反褶积去噪结果的频谱相比频带更宽，对高低频信号能量保护性更好，并且与原始信号的频带也更为接近.
|
图 7 两种方法去噪结果频谱分析对比 Fig. 7 Comparison of two methods for denoising results spectrum analysis |
本文提出利用Hessian矩阵特征值对应线性目标的关系在多个尺度上对随机噪声进行压制的方法，通过理论模型测试和实际资料处理，可以得到以下结论和认识.
(1) 基于Hessain矩阵特征值的多尺度去噪算法可以利用有效信号与随机噪声对应Hessian矩阵的特征值差异检测出剖面中的线性目标，从而将有效信号与随机噪声分离开来，达到去噪效果.
(2) 小尺度对应高频信号，大尺度对应低频信号，要想兼顾有效信号低频和高频的能量必须在多个尺度上联合进行去噪处理.
(3) 通过与传统的F-X反褶积去噪算法对比可见，本文方法去噪效果更好，并能保持有效波的完整性.
线性滤波器的响应函数是影响本文方法去噪结果的关键，寻找性能更好、对有效信号更敏感的滤波器的响应函数是我们需要进一步研究的内容.
Gulunay N. 1986. FXDECON and complex wiener prediction filter.//56th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 279-281.
|
Yu S, Cai X, Su Y. 1989. Seismic signal enhancement by polynomial fitting. Applied Geophysics, 1(1): 57-65. |
Zhang J H, Lv N, Tian L Y, et al. 2005. Comprehensive review of seismic data denoising. Oil Geophysical Prospecting (in Chinese), 40(S1): 121-127. |
Liu W, Yang K, Wang Z, et al. 2010. Orthogonal polynomials transform based on SVD and its application in seismic dataprocessing. Chinese Journal of Engineering Geophysics (in Chinese), 7(4): 433-437. |
Wang Y, Li J, Zhang K F. 1998. The noise elimination processing of seismic signal based on wavelet transform. Geophysical Processing for Petroleum (in Chinese), 37(3): 72-76. |
Frangi A F, Niessen W J, Vincken K L, et al. 2000. Multiscale vessel enhancement filtering.//Wells W M, Colchester A, Delp S eds. Medical Image Computing and Computer-Assisted Intervention. Berlin, Heidelberg: Springer, 1496: 130-137.
|
Koller T M, Gerig G, Szekely G, et al. 1995. Multiscale detection of curvilinear structures in 2-D and 3-D image data.//Proceedings of IEEE International Conference on Computer Vision. Cambridge, MA, USA, USA: IEEE, 864-869.
|
Lorenz C, Carlsen I C, Buzug T M, et al. 1997. Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2D and 3D medical images.//Troccaz J, Grimson E, Mösges R eds. CVRMed-MRCAS'97. CVRMed 1997, MRCAS 1997. Berlin, Heidelberg: Springer, 1205: 233-242.
|
Xu X, Liu B, Zhou F G. 2013. Hessian-based vessel enhancement combined with directional filter banks and vessel similarity.//2013 ICME International Conference on Complex Medical Engineering. Beijing, China: IEEE, 80-84.
|
Shashank, Bhattacharya M, Sharma G K. 2013. Optimized coronary artery segmentation using Frangi filter and anisotropic diffusion filtering.//2013 International Symposium on Computational and Business Intelligence. New Delhi, India: IEEE, 261-264.
|
You J, Chen B. 2011. New approach to retinal image enhancementbased on Hessian matrix. Journal of Computer Applications, 31(6): 1560-1562. DOI:10.3724/SP.J.1087.2011.01560 |
Qin H, Huang X X. 2016. Coronary angiography image segmentation and skeleton extraction based on Hessian matrix. Journal of Data Acquisition & Processing, 31(5): 911-918. DOI:10.16337/j.1004-9037.2016.05.007 |
Koenderink J J. 1984. The structure of images. Biological Cybernetics, 50(5): 363-370. DOI:10.1007/BF00336961 |
Florack L M J, Romeny B M T H, Koenderink J J, et al. 1992. Scale and the differential structure of images. Image and Vision Computing, 10(6): 376-388. DOI:10.1016/0262-8856(92)90024-W |
张军华, 吕宁, 田连玉, 等. 2005. 地震资料去噪方法综合评述. 石油地球物理勘探, 40(S1): 121-127. DOI:10.3321/j.issn:1000-7210.2005.z1.025 |
刘伟, 杨凯, 王征, 等. 2010. 基于SVD的正交多项式变换及其在地震资料处理中的应用. 工程地球物理学报, 07(4): 433-437. DOI:10.3969/j.issn.1672-7940.2010.04007 |
王勇, 张奎凤, 郦军. 1998. 基于小波变换的地震信号降噪处理. 石油物探, 37(3): 72-76. |
游嘉, 陈波. 2011. 基于Hessian矩阵的多尺度视网膜图像增强方法. 计算机应用, 31(6): 1560-1562. DOI:10.3724/SP.J.1087.2011.01560 |
秦红星, 黄晓雪. 2016. 基于Hessian矩阵的冠脉造影图像分割与骨架提取. 数据采集与处理, 31(5): 911-918. DOI:10.16337/j.1004-9037.2016.05.007 |