0 引言

众所周知，地震勘探是勘探地下油气等资源的重要方法.地震勘探资料通常会伴随各种各样的随机噪声，这些噪声势必会对勘探中各种重要数据的采集产生一定的干扰.其中，激发前就存在的背景噪声(即初至前噪声)是不可避免的(徐义贤和罗银河，2015；田雅男等，2015).因此，为了能够有效提取数据，就要尽可能去了解地震勘探过程中的各种随机噪声从而对噪声进行压制(林红波等，2011).在实际地震勘探中，环境噪声占有很大比重，而风成环境噪声又是环境噪声的主要成分(彭蛟龙，2013；李光辉和李月，2016).因此，有必要对风成环境系统进行建模研究.

目前，关于地震勘探中环境噪声的研究大多是对其性质进行探讨，而对环境噪声的成因研究甚少，并且几乎没有关于地震勘探环境系统建模的研究与分析.本文对沙漠地区地震勘探风成环境系统的建模研究，其结果正是对环境噪声的成因说明，也是对地震勘探环境系统建模的首探研究.风成环境噪声的主要作用源是风.风作用在检波器上的微震(Astiz and Creager, 1994; Chávez-García and Rodríguez, 2007)以及在勘探地区的风吹草动等等都会以噪声的形式反映在所采集的数据中，形成干扰(Li et al., 2015).1996年，Mitchell M.Withers等人提出风速与地震背景噪声之间具有很强的相关性(Wither et al., 1996)，并进行了研究.2005年，Marco Mucciarelli等人分析了风对于微震测量的影响，认为风会增加微震波场中所有分量的振幅，并且空气，仪器和地面接口处中的气流会对实验中的传感器或数字检波器产生不利的影响(Mucciarelli et al., 2005).2011年，Barajas-Olade和Ramadan利用现场测量的方法研究了随机噪声与风速之间的相关性(Barajas-Olade and Ramadan., 2011).以上均可说明，风与地震勘探中的环境噪声存在一定的关系.

本文主要从以下三个方面进行研究.首先，采用自回归滑动平均(ARMA)模型对随机风序列进行模拟，模拟得到的结果视为风成环境系统的输入；然后，对沙漠地区风成环境系统的性质进行分析，得知勘探地区所耦合成的系统具有混沌特性；首次采用duffing方程对风成环境系统进行建模，将模拟得到的随机风序列输入到系统中，将输出信号与沙漠地区初至前噪声记录分别从时域波形，相态图以及频谱图三个方面进行对比.结果显示，我们所建立的沙漠地区风成环境系统是可行且有效的.

1 自回归滑动平均(ARMA)模型模拟随机风序列

由于沙漠地区所对应勘探时刻的风速数据未知，因此需要对风速序列进行模拟.风速序列是一种随机事件序列，其变化的随机性很强.而对于这种随机时间序列的分析，自回归滑动平均模型在建模和预测方面具有较好的效果.根据风致振动理论，对风速序列的模拟主要是完成对其成分中脉动风的模拟，再叠加平均风速即可(Li et al., 2015).较低阶的ARMA模型，(一般在(4, 4)阶以下)(李英民等, 1993, 1992)即可较好地模拟风序列中的脉动风.确定模型阶数后，利用自相关函数方法估计模型参数，得到最终风速序列.

1.1 ARMA模型

对于一个具有加权形式的ARMA(n, m)模型的一般形式为

(1)

其中，{x_t}是随机时间序列，均值为零；φ_n≠0, 是自回归系数，n是自回归阶数；θ_m≠0，是滑动平均系数，m是滑动平均阶数；{ε_t}是白噪声过程.

自回归滑动平均ARMA(n, m)模型的自相关函数可以看成是滑动平均MA(m)模型的自相关函数和自回归AR(n)模型的自相关函数的混合.经过推导，当k>m时，ARMA(n, m)的自相关函数为

(2)

其中，ρ_k由观测数据得到(在本文中，ρ_k是Davenport谱的傅里叶逆变换).根据式(2)，通过一组列方程组，即可求得自回归系数.然后改写模型，令式(1)中，x′_t=x_t－φ₁x_t－1－φ₂x_t－2－…－φ_nx_t－n，则式(1)变为

(3)

此式实际上为MA(m)模型.

推导得到式(3)的自相关函数为

(4)

式中，ρ′_k是x′_t的自相关函数(Wei, 1994; Adas, 1997; Moon and Stirling, 2000; 徐文令等，2013).根据式(4)，通过一组列方程组，即可求得滑动平均系数.

1.2 达文波特(Davenport)风谱

脉动风速谱描述的是紊流风的谱特性，用以分析频谱内的随机响应，是用来对脉动风特性进行描述的一个十分重要的特征.本文选取了被广大研究学者高度认可的脉动风速谱中的达文波特(Davenport)谱，用以拟合模拟出来的随机风序列的功率谱，以验证所模拟的精确程度.经验公式为

(5)

式中, S_v(f)表示脉动风速谱，f为脉动风频率；K为地面粗糙系数；v₁₀为距地面10 m高度处的平均风速；u为湍流积分尺度系数(张立超，2014；Li et al., 2015).当自回归阶数n为2，滑动平均阶数m为1，K=0.005，v₁₀=16 m/s时，求得的随机风序列(见图 1)的功率谱与Davenport谱的拟合如图 2所示.如图可知，拟合程度较好，证明利用ARMA模型可以有效模拟随机风序列.

2 风成环境系统建模

系统建模就是根据某个特定目的通过一定的数学方法对某系统或其子系统进行建模(王涛，2008)，并建立系统数学模型.数学模型就是描述系统内、外部变量间相互关系的数学表达式，它的最终目的是确定系统模型形式、结构和参数，获得正确反映系统表征、特征和功能的数学表达形式(陈腾浩等，2012).本文将沙漠勘探地区的地形、地表和地质等条件以及勘探位置与检波器的耦合程度作为影响因子，构成沙漠地区风成环境系统，并对其进行建模研究.

2.1 脉冲响应法分析系统性质

勘探地区风成环境系统比较复杂，为准确建立模型，须先简单分析系统的性质.定性分析，风速序列是具有混沌性质的随机序列(Hara et al., 1997)，而初至前噪声也是具有混沌性质的序列(彭蛟龙，2013).因此，可初步判定:该系统可以是线性系统，也可以是混沌系统.

因其系统较为复杂，为方便分析，先采用非参数模型法(侯忠生，1999；Ferraty and Vieu, 2006)进行分析.非参数模型是指从系统的实验过程中，直接或间接所获得的响应，可以利用记录或分析系统的输入和输出信号来估计系统.假设勘探地区风成环境系统为线性定常系统，采用脉冲响应法来求得风成环境系统的脉冲响应.脉冲响应法(方崇智和萧德云，1998)是利用被辨识线性定常系统的输入和输出信息，通过脉冲响应来辨识系统的数学模型.

已知线性定常系统卷积积分式为

(6)

其中，u(t)为该系统输入；y(t)为系统输出；g(t)为该系统的脉冲响应；假定u(t)和y(t)被一个采样周期为T的采样器周期性采样，当T足够小时，逐段常值逼近可将u(t)、y(t)和g(t)用阶梯信号近似代替，公式为

(7)

推导可得：

(8)

其中，

因此，式(8)可以表示为

(9)

综上可知，系统的脉冲响应为

(10)

实验中，将由ARMA(2, 1)模型模拟得到的随机风序列作为输入，提取沙漠某勘探区第1道1000点初至前噪声记录作为输出，如图 3a.根据上述方法，求得系统的相态图，如图 3b.

混沌序列的相态图表现为在有限空间内不断伸长和折叠的回复性运动(彭蛟龙，2013).从相态图 3b可知，该系统具有混沌性质，即实验结果与预先假设的线性定常性质不符，因此利用混沌特性对风成环境系统进行建模分析，本文将采用具有混沌性质的duffing方程模型对系统进行建模.另外，从地球本身来说，地球是一个耗散系统，耗散系统中的混沌吸引子说明地球系统是一个非线性系统(刘式达，1990)；而地震勘探环境系统可视为缩小的地球系统，这从侧面间接的验证了我们的证明.

2.2 duffing方程模型模拟沙漠地区风成环境系统

终上所述，地震勘探风成环境系统可视为是具有混沌特性的系统.由此，选用如下duffing方程作为系统模型，公式为

(11)

考虑到沙漠地区初至前噪声具有类余弦周期的波动特性(如图 3a所示), 因此在方程右侧的驱动项中填加余弦控制项.其中，k₁为阻尼比(Albarello and Baliva, 2009)；k₂为线性回复力系数；k₃和k₄为非线性回复力系数；r为策动力幅度；ω为策动力频率；t为采样点时间；st为系统输入，即随机风序列.r和ω控制输出信号y的幅度和类周期状况；当方程右侧给出不同的策动力时，方程输出就会呈现不同的状态(李月等，2005; Sheu et al., 2007; Ma and Zhang, 2010；Kovacic and Brennan, 2011；刘财等，2015).因此，根据不同参数对duffing方程输出结果的影响，将得到的duffing方程输出信号与勘探地区初至前噪声记录进行对比，就可确定duffing方程的各个参数.

2.3 duffing方程参数分析

在实际地震勘探中，勘探地区不同道的初至前噪声的时域波形不可能是完全相同的，但其频谱和相态图的差别不大，其走向和规律是相同的(彭蛟龙，2013).重要的是，频谱和相态图表现了环境噪声的频率成分和混沌特性(Campillo, 2006; 彭蛟龙，2013).因此，在分析方程参数对输出结果的影响时，将从道记录的时域波形，相态图和频谱图三个方面来进行对比.由于无法得知勘探地区采集时刻初至前噪声记录所对应风序列的风速大小，所以无法保证依据davenport谱所模拟出来的随机风序列的风速大小与实际相符，因此，我们将分别采用归一化后的时域波形，相态图和频谱图来进行对比.

要保证duffing方程输出具有混沌特性，各个参数具有一定的选取范围，并且各个参数相互制约.在实验时，对于duffing方程中的六个参数，为验证一般性，令r=ω=1，剩余四个参数，根据经验参数设定三个，变化另外一个，观测时域波形，相态图和频谱图的变化，分别记录四个参数的变化规律，方便后续参数调整.然后对r进行调整，过程同上.最后，通过调整ω来控制duffing方程输出结果的类周期宽度.

记录结果如下：在除ω外的五个参数中，k₁和r对duffing方程输出结果影响较大，控制着输出信号时域波形上的毛刺状况，相态图的形状和覆盖范围以及频谱幅度大小.随着k₁的增大，输出信号时域信号波形上的毛刺减少从而波动形态逐渐清晰；相态图的形状和覆盖范围发生变化.当k₁增大到一定程度时，这些变化不再明显(参考示例，如图 4).与k₁类似，随着r的增大，输出信号时域信号波形上的毛刺减少从而波动形态逐渐清晰；相态图的形状和覆盖范围发生变化，除此之外，输出信号的频谱幅度增大.但当r增大到一定程度时，输出信号时域波形上的毛刺增多从而波动形态模糊，且频谱幅度降低(参考示例，如图 5).k₂，k₃和k₄对输出结果影响较小，但当线性回复力k₂增大到一定程度时或者当非线性回复力k₃和k₄同时为0时，duffing方程输出结果将不再混沌.另外，在选择参数的过程中，还要兼顾不同参数之间的相互影响，改变一个参数的同时对其他参数进行微调以达到预期结果.

以下给出对duffing方程输出结果影响较大的参数k₁和r的调节过程示例(实际参数调节中，工作量较大，因此本文仅给出代表性示例)：

(1) 固定k₂=2.5，k₃=15，k₄ =3，r=ω=1，k₁分别取0.01，10，30.结果依次如图 4所示.

(2) 固定k₁=23，k₂=5，k₃=20，k₄=0.5，ω=1，r分别选取-1，-0.01，1，2，5.结果依次如图 5所示.

选取上述1000点沙漠地区初至前噪声记录，其时域波形，相态图和频谱如图 6a、c和e所示.仿真过程中，通过大量实验，多次调整，并不断与图 6a、c和e进行拟合对比，可得到当参数k₁=27，k₂=25，k₃=20，k₄=0.55，r=2.6时，与沙漠地区初至前噪声记录拟合程度最好，然后通过控制ω以调整duffing方程输出信号的类周期宽度尽可能与沙漠地区实际记录吻合，经多次实验，ω=0.05时效果最佳.图 6b、d和f为duffing方程输出结果.

如图 6所示，从整体上看，duffing方程输出信号的时域波形上的毛刺状况，波动形态和类周期宽度均与沙漠地区初至前噪声记录相吻合；相态图重叠范围较吻合，但整体趋势有所偏差，实际记录整体向左微偏，而duffing方程输出信号整体趋势垂直于横轴，此处有待继续改进；输出信号频谱幅度稍高于初至前噪声幅度，但频谱截止频率均在50 Hz左右，十分吻合.以上证明了采用混沌性duffing方程模拟风成环境系统具有一定的可行性和有效性.

3 结论

本文首先提出了ARMA模型方法模拟随机风序列.然后，分析得知，沙漠地区风成环境系统具有混沌性质，并采用duffing方程模型模拟沙漠地区风成环境系统.以随机风序列作为系统激励，通过系统后的输出信号在时域波形，相态图和频谱三个方面与沙漠地区初至前噪声记录结果拟合较符，证明了采用duffing方程模拟风成环境系统具有一定的可行性和有效性.结果表明，初至前噪声记录是由风作用于勘探地区所形成的，即风是产生实际陆地地震勘探中环境噪声的主要原因，对现有环境噪声的认识进行了补充，为以后进一步提高地震勘探去噪技术提供参考依据，从而获得更高质量的实际地震记录，精确找到地下石油和天然气等资源.

致谢 感谢国家自然科学基金(41130421和41574096)的资助，感谢编辑和审稿人对本文的编辑与审阅.