基于噪声建模分析的沙漠地震勘探随机噪声消减方法与流程

本发明属于地球物理技术领域，针对传播介质的非均匀性对波场传播物理特征的复杂表现，分析随机噪声成因并建立随机噪声模型，最终利用随机噪声模型来指导选择合适的噪声消减方法。

背景技术：

在地震勘探信号处理中，随机噪声严重影响了有效信号的获取，因此对于随机噪声的研究和建模是信号处理的首要步骤。通过对随机噪声性质分析及噪声建模，可以调整模型参数模拟不同条件下的随机噪声进而针对不同的噪声选取合适的消噪方法。随机噪声模型越接近实际噪声，对于噪声性质的分析和滤波方法选取就越准。finette和heimeyer(1996)建立了风关噪声场模型，并指出必须考虑更复杂的环境参数如流体静压，海浪与气流，涡流，空气边界层的响应等。在最近十年对陆地随机噪声的研究中，我们知道短时长陆地地震勘探随机噪声不满足高斯随机噪声，也不能简单地归为线性随机噪声且受到采集环境的影响很大。相比于常规噪声，地震勘探随机噪声性质更加复杂以至于常用的滤波方法无法达到令人满意的去噪效果，亟需进一步研究为实际消噪处理提供可靠的理论基础。目前已有的建模技术，如li(2017)等人提出了在均匀介质中的随机噪声建模，其中假设噪声传播介质是均匀各向同性的理想介质。这个模型是学者对于随机噪声建模的初步尝试，其假设条件过于理想与实际条件有较大差距，忽略了实际传播介质中物理因素对传播特征带来的复杂表现，会对后续去噪处理带来较大误差，因此存在很大的改进空间。

技术实现要素：

本发明提供一种基于噪声建模分析的沙漠地震勘探随机噪声消减方法，以解决现有的建模方法忽略了实际传播介质的复杂物理因素从而不能全面描述传播特征的问题，从而为后续消噪方法的选择提供有效的依据和必要的合成数据。针对非均匀介质对传播速度的改变，建立与实际噪声高度拟合的随机噪声模型，为后续信号处理和噪声压制提供全面的理论基础，最终利用随机噪声模型构成的合成记录指导选择适合于沙漠地震勘探随机噪声的消噪方法。

本发明采取的技术方案是，包括下列步骤：

(1)构建含噪合成记录，沙漠地震勘探合成记录为检波器记录的地下反射信号即有效地震勘探信号和环境噪声等随机噪声的叠加，表示为：

g(d,t)＝s(d,t)+φ(d,t)(1)

其中d和t分别表示道号和时间，g(d,t)表示含噪信号，s(d,t)表示有效地震勘探信号，s(d,t)地震子波采用雷克子波构成，雷克子波定义如下：

s(t)＝(1-2π²f²t²)exp(-π²f²t²)(2)

其中f为中心频率；

φ(d,t)为沙漠地震勘探随机噪声模型；

(2)构建沙漠地震勘探随机噪声模型φ(d,t)

测区为沙漠陆地地区，将随机噪声分为三类：自然噪声，近场人文噪声，远场人文噪声，并假设噪声源以点源的形式随机分布在检波器周围；

在沙漠地区植被稀少且地面较平坦，自然噪声主要是由于风吹地面产生的，则根据风阵理论的风压表达式，得到自然噪声源的表达式为：

其中是风压谱的第i个分量，ρ是空气密度，为高度为sm时的风速，s(w)为脉动风速谱w角频率，k是地表拖曳力系数，为参考风速，χ为湍流积分尺度系数，且满足wi是第i个分量的角频率，εi是第i个分量的初始相位，是一个均匀分布在(0，2π)的随机变量，n是分量个数，是角频率增量，wh和wl分别是脉动风谱能量分布频率范围的上限和下限，t是时间变量；

近场人文噪声主要由机器工作产生，利用伪高斯函数模拟近场噪声源函数：

其中a为振幅，bw为带宽，τ为时间延迟，fc为中心频率，为随机初始相位；

远场人文噪声主要是由检波器远处的公路，工厂和城镇引起的噪声，远场人文噪声源其表达式如下(lietal.2017)：

其中aj是振幅，kj是波数，wj是角频率，为随机初始相位，x为位移；

此时波在弱非均匀介质中传播，设弱非均匀介质中拉梅系数λ＝λ0(1-0.001x)，μ＝0，密度ρ为常数；

波动方程描述了波场传播规律，非均匀介质中的波动方程：

其中其中为波场函数，f为激励函数，c(r)是非均匀介质中传播速度，r是波场传播位移；

将三个噪声源函数分别作为波动函数的激励函数，即式(3)，(4)，(5)分别作为式(6)的右端项f(r,t)，在初始条件和边界条件下，用格林函数求解波动方程，分别得到自然噪声的噪声场为:

其中nw为噪声点源个数，n是风压谱的分量个数，r为点源位置矢量，rj为点源的位置向量，cn是噪声在非均匀介质的传播速度，cn＝c0(1-0.001x)，εi是在0～2π内随机变化的初始相位；

近场人文噪声场为：

其中nw为噪声点源个数，vn为非均匀介质中纵波传播速度vn＝v0(1-0.001x)；

远场人文噪声场为：

其中nx为检波器的个数，vk为非均匀介质中噪声传播速度vk＝v0(1-0.001x)；

则三种噪声场的叠加得到沙漠地区随机噪声模型的表达式为：

φ(r,t)＝φw(r,t)+φm(r,t)+φf(r,t)(10)

φw(r,t)，φm(r,t)，φf(r,t)分别是式(7)自然噪声场，式(8)近场人文噪声场和式(9)远场人文噪声场的表达式；

这样，模拟d道的噪声模型记录，得到φ(d,t)，得到含噪合成记录即式(1)；

(3)利用随机噪声模型确定适合的滤波方法进行滤波消噪

采用shearlet消噪方法压制沙漠地震勘探随机噪声，shearlet变换阈值消噪的步骤如下：

shearlet变换是线性的，因此式(1)中含噪记录g(d,t)被变换至shearlet域，其shearlet域系数定义为：

shg(d,t)＝shs(d,t)+shφ(d,t)(11)

其中sh是shearlet变换系数，然后对含噪信号进行shearlet变换阈值去噪，得到滤波消噪后的信号为：

其中，t是阈值算子并且为硬阈值：其中a表示信号，是消噪后的信号，其中sh和sh^-1分别是shearlet变换和反变换系数。

本发明的优点：建立了与实际噪声更加相似的随机噪声模型并借助随机噪声模型构成的含噪记录调整消噪程序中的关键参数，进而在分析噪声性质和后续选择有效消噪方法中达到更加准确可靠的效果，最终根据随机噪声模型提供的先验知识和合成数据的实验结果确定了适合于沙漠地震勘探随机噪声的消减方法即shearlet变换法。本方法在考虑到非均匀介质给传播性质带来的复杂表现，使随机噪声模型与实际噪声的数理特征更为相似，达到更好的拟合效果。而准确的随机噪声模型和分析是信号处理的基础，为后续信号处理提供了准确的先验知识和可靠的数理基础，通过随机噪声模型选择合适的滤波方法并确定重要参数进而推广到实际信号处理中，本发明基于随机噪声模型选择了shearlet消噪方法，shearlet变换后能量的分布会随着变换尺度的变换尺度呈现出一定的规律，会出现相对的低频部分和高频部分，低频部分包含了原始图像的大部分能量，是原图像在不同尺度下的逼近信号。通过对基本函数得到膨胀、剪切和平移shearlet变换具有良好的尺度特性、方向敏感性以及局部化特性，对于沙漠地震随机噪声与信号混叠，干扰强度大等特点有良好的去噪效果。

附图说明

图1(a)是由风引起的自然噪声源的分布区域图，风噪声源以点源的形式随机且不均匀地分布在检波器左侧三角形区域中，根据风源影响范围将风源作用区域边长设为500m；

图1(b)是近场人文噪声源的分布区域图，根据实际工作经验，设定近场人文的噪声源随机分布在半径为200m到1250m的圆形区域内。

图1(c)是远场人文噪声源的分布区域图，远场人文噪声源随机且不均匀分布在距检波器1500m处，由于距检波器距离较远可以看做距离是一定的因此噪声源呈线状分布；

图2(a)是在实际地震资料中截取的一块沙漠地区地震勘探随机噪声图；

图2(b)是利用本发明建模方法根据图2(a)实际环境调整模型参数建立的随机噪声模型；

图2(c)在均匀介质中建立的随机噪声模型，可以看出均匀介质中随机噪声模型的波型变化比较简单，而本发明模拟的沙漠地震勘探随机噪声在波形上与实际噪声极其相似；

图3(a)是图2(a)的fk谱图，用来表示信号的能量分布；

图3(b)是图2(b)的fk谱图；

图3(c)是2(c)的fk谱图；可以看出均匀介质中的随机噪声模型能量分布范围与实际噪声存在差别，比实际噪声频率要低一些，本发明建立随机噪声模型与实际随机噪声模型在能量分布上基本一致；

图4(a)是对图2(a)和图2(b)中50道数据分别求功率谱密度然后再取其平均得到的平均功率谱密度，实线代表本发明建立随机噪声模型图2(b)，虚线代表实际噪声记录图2(a)；

图4(b)是对图2(c)和图2(a)求平均功率谱密度，实线代表均匀介质中建立的随机噪声模型图2(b)，虚线代表实际噪声记录图2(a)，可以看出本发明建立随机噪声模型和实际随机噪声模型在能量强度，分布范围和主要能量集中范围都非常相似，均匀介质中随机噪声模型的频率分布范围与实际的存在差别；

图5(a)是分别在图2(a)和(b)随机取一道数据(第45道)计算得到累积概率分布图，实线代表本发明建立随机噪声模型图2(b)，虚线代表实际噪声记录图2(a)；

图5(b)是分别在图2(a)和(c)随机取一道数据(第45道)计算得到累积概率分布图，实线代表均匀介质中建立的随机噪声模型图2(b)，虚线代表实际噪声记录图2(a)；计算结果完整描述了两个随机噪声的概率分布，可以看出本发明建立随机噪声模型与实际噪声的概率分布曲线几乎吻合；

图6(a)是含噪合成地震记录，其中噪声是本文建立的随机噪声模型；

图6(b)shearlet消噪结果；

图6(c)curvelet消噪结果，从两个消噪结果可以看出，curvelet消噪方法对随机噪声压制效果较差；shearlet消噪后，将随机噪声和面波都得到有效压制，同相轴恢复清晰连续，说明对于模拟的随机噪声模型得到的合成记录，shearlet方法对能够有效压制沙漠地震勘探随机噪声；

图7(a)是实际沙漠地震资料；

图7(b)shearlet消噪结果；

图7(c)curvelet消噪结果，从消噪结果上可以看出，对于实际地震资料，curvelet消噪去除少部分噪声，大部分随机噪声和面波都保留着，消噪效果较差，而shearlet消噪后，大量随机噪声和面波都被压制，有效信号清晰连贯，可以看出对于实际地震资料消噪效果与合成资料一致，说明通过随机噪声模型选择的滤波方法是准确的，这对我们选择其他的滤波方法起到了很好指导作用。

具体实施方式

包括下列步骤：

(1)构建含噪合成记录，沙漠地震勘探合成记录为检波器记录的地下反射信号即有效地震勘探信号和环境噪声等随机噪声的叠加，表示为：

g(d,t)＝s(d,t)+φ(d,t)(1)

其中d和t分别表示道号和时间，g(d,t)表示含噪信号，s(d,t)表示有效地震勘探信号，s(d,t)地震子波采用雷克子波构成，雷克子波定义如下：

s(t)＝(1-2π²f²t²)exp(-π²f²t²)(2)

其中f为中心频率；

φ(d,t)为沙漠地震勘探随机噪声模型；

(2)构建沙漠地震勘探随机噪声模型φ

测区为沙漠陆地地区，将随机噪声分为三类：自然噪声，近场人文噪声，远场人文噪声，并假设噪声源以点源的形式随机分布在检波器周围，如附图1所示；

在沙漠地区植被稀少且地面较平坦，自然噪声主要是由于风吹地面产生的，则根据风阵理论的风压表达式，得到自然噪声源的表达式为：

其中是风压谱的第i个分量，ρ是空气密度，为高度为sm时的风速，s(w)为脉动风速谱w角频率，k是地表拖曳力系数，为参考风速，χ为湍流积分尺度系数，且满足wi是第i个分量的角频率，εi是第i个分量的初始相位，是一个均匀分布在(0，2π)的随机变量，n是分量个数，是角频率增量，wh和wl分别是脉动风谱能量分布频率范围的上限和下限，t是时间变量；

近场人文噪声主要由机器工作产生，利用伪高斯函数模拟近场噪声源函数：

其中a为振幅，bw为带宽，τ为时间延迟，fc为中心频率，为随机初始相位；

远场人文噪声主要是由检波器远处的公路，工厂和城镇引起的噪声，远场人文噪声源其表达式如下(lietal.2017)：

其中aj是振幅，kj是波数，wj是角频率，为随机初始相位，x为位移；

此时波在弱非均匀介质中传播，设弱非均匀介质中拉梅系数λ＝λ0(1-0.001x)，μ＝0，密度ρ为常数；

波动方程描述了波场传播规律，非均匀介质中的波动方程：

其中其中为波场函数，f为激励函数，c(r)是非均匀介质中传播速度，r是波场传播位移；

将三个噪声源函数分别作为波动函数的激励函数，即式(3)，(4)，(5)分别作为式(6)的右端项f(r,t)，在初始条件和边界条件下，用格林函数求解波动方程，分别得到自然噪声的噪声场为:

其中nw为噪声点源个数，n是风压谱的分量个数，r为点源位置矢量，rj为点源的位置向量，cn是噪声在非均匀介质的传播速度，cn＝c0(1-0.001x)，εi是在0～2π内随机变化的初始相位；

近场人文噪声场为：

其中nw为噪声点源个数，vn为非均匀介质中纵波传播速度vn＝v0(1-0.001x)；

远场人文噪声场为：

其中nx为检波器的个数，vk为非均匀介质中噪声传播速度vk＝v0(1-0.001x)；

则三种噪声场的叠加得到沙漠地区随机噪声模型的表达式为：

φ(r,t)＝φw(r,t)+φm(r,t)+φf(r,t)(10)

φw(r,t)，φm(r,t)，φf(r,t)分别是式(7)自然噪声场，式(8)近场人文噪声场和式(9)远场人文噪声场的表达式；

这样，模拟d道的噪声模型记录，得到φ(d,t)，得到含噪合成记录即式(1)；

(3)利用随机噪声模型确定适合的滤波方法进行滤波消噪

通过建立必要的合成含噪记录(1)进行实验才能选择有效的滤波方法，针对沙漠随机噪声与信号频带混叠不易去除的问题，我们选择换域去噪方法shearlet并利用随机噪声模型构建的合成含噪记录进行实验，最终确定shearlet消噪方法更适合压制沙漠地震勘探随机噪声，下面是shearlet变换阈值消噪方法的步骤：

shearlet变换是线性的，因此式(1)中含噪记录g(d,t)可以被变换至shearlet域，其shearlet域系数定义为：

shg(d,t)＝shs(d,t)+shφ(d,t)(11)

其中sh是shearlet变换系数，然后对含噪信号进行shearlet变换阈值去噪，得到滤波消噪后的信号为：

其中，t是阈值算子并且为硬阈值：其中a表示信号，是消噪后的信号，其中sh和sh^-1分别是shearlet变换和反变换系数；

通过含噪合成记录g(d,t)调整阈值算子t并得到最优阈值t1，使shearlet变化消噪方法对合成含噪记录的去噪效果达到最优，进而这个将阈值t1用到实际地震资料的消噪中达到良好的去噪效果。

下边通过具体实验来进一步说明本发明。

实验例1本发明建立随机噪声模型的准确性

利用本发明提出随机噪声模型建立一个50道的模拟噪声记录，每道800个采样点，道间距15m。如图2(b)所示，并与图2(a)真实沙漠随机噪声记录对比，可以看出本发明模拟随机噪声模型和实际噪声波型基本一致，并且相比于在均匀介质中建立的随机噪声模型图2(c),本发明提出噪声与实际噪声更加相似。图3是对图2中3幅记录分别求fk谱图得到的结果，可以看出在能量分布上，本发明提出随机噪声模型与实际基本一致。图4是对图2中每一幅记录求平均功率谱密度，可以看出本发明提出随机噪声模型在能量大小和主频分布范围上都与实际噪声基本一致同时也发现沙漠地震勘探随机噪声的频率很低，其主频分布大约在10hz左右。图5是随机抽取一道记录(本发明是第45道)，并计算其累计概率分布，从对比结果可以看出，本发明随机噪声模型的累计概率分布曲线与实际基本吻合，而在均匀介质中的随机噪声模型的概率分布与实际噪声有差别。综上所述本发明提出随机噪声模型与实际噪声更加一致。以上对比实验从时间域和空间域验证了本发明提出随机噪声模型与实际噪声数理特征基本吻合，证明了本文建模方法的合理性。

实验例2利用随机噪声模型确定适合的滤波方法

为了证明本发明构成的含噪合成记录可以用于选择滤波算法，选取两种算法应用于合成地震勘探记录中。合成记录中包含两个不同程度的弯轴和两个断轴。道间距为15m,采用频率为800hz。浅层地震子波主频为25hz，速度为1800m/s,深层地震子波主频为20hz,速度为2000m/s。两个断轴地震子波主频为20hz,速度为1800m/s。将本发明建立的随机噪声模型加入到合成记录中得到含噪合成记录如附图6(a)所示。

选取去噪方法shearlet应用于通过随机噪声模型合成的含噪记录,并将curvelet作为对比实验。通过消噪效果对比来选择适合于沙漠地震勘探随机噪声的消噪方法。

结果分析:从附图6中可以看出对于含噪合成记录，shearlet对于沙漠地震勘探随机噪声有着良好去噪效果，大量随机噪声被压制且同相轴恢复清晰，而curvelet对于此类低频随机噪声基本没有消噪效果。进而推广到对实际地震资料的处理如附图7所示，可以看出消噪结果对比与合成记录处理结果一致，shearlet对于沙漠地震勘探随机噪声的消噪效果明显优于curvelet。最终确定shearlet变换阈值去噪方法适合于沙漠地震随机噪声压制工作。

综上所述，消噪结果的对比体现了由随机噪声模型构成的含噪合成记录对于选择去噪方法的指导作用以及在消噪实验中的应用。通过对于噪声的建模和分析，我们可以了解噪声源一系列数理特征并且构成实验所必须的含噪合成记录用于消噪实验，从而针对性的选择合适的滤波方法，本发明基于噪声建模最终确定shearlet变化阈值去噪方法适合于沙漠地震勘探噪声的消减，能够提高地震勘探资料信噪比，更准确地恢复有效信号，有效地帮助了后续地震资料解释。