随着地球物理勘探技术的发展，具有安全可控、绿色环保等特点的各类非爆炸震源得到了广泛的应用^[1]。以液压式震源为代表的大型人工震源设备在地表平坦、远离人烟的地区已取得较好的勘探效果^[2]，但在人口密集、森林覆盖率较高的区域以及丘陵山地，急需一种结构简单、运移方便且有较大输出力的小型震源装置来弥补这一空白。精密可控震源是20世纪末发展起来的一种新型震源，其结构简单、易于拆装，噪音小且无尾气排放，不仅符合目前各国绿色探勘对震源设备的要求^[3]，同时也十分适合在大型震源无法进入的特殊地区进行勘探作业。缘于精密可控震源本身结构特点，它激发的扫描信号具有重复性好、频率低、传播距离远等优点，为石油、天然气探勘开发提供了一条新路径^[4-7]。

震源装置信号激发是地震勘探的第一个环节，为提升地质反演资料的分辨率，以期最大程度地将振动器产生的相干振动能量传入地下，且尽量减少输出信号的畸变^[8]，从而更加准确地推断地下油气资源分布规律。由于精密可控震源压重通过隔振系统与平板连接在一起，在实际激振作业过程中，激振器带动平板以及大地表层振动的同时不可避免地使压重产生振动。由于压重本身的质量较大，它产生的振动不仅会削减输出信号的能量，同时会产生干扰噪声，降低震源输出信号质量。不合理的隔振系统参数会导致压重振幅过大，甚至产生共振。

目前针对动力机械上部隔振系统已有了不少研究：2001年牛军川等基于能量传递观点分析了橡胶隔振器对旋转机械的隔振效果^[9]；2009年邓东等针对振动器阀控系统建立了动力学模型并求解，寻找最优阻尼比以改善扫描起始段信号质量^[10]；2013年李兆锋等采用有限元方法分析了双层浮筏隔振系统在振动筛中的隔振效果^[11]；2014年李淑清等针对液压式可控震源低频段设计了混合隔振系统^[12]；2015年王永芳等通过有限元分析了隔振参数对液压式可控震源上装部件的影响规律^[13]。但是这些研究往往忽略了大地对上部结构振动特性的影响，所求得的结论可能会与实际情况有所偏差。因此，在分析隔振系统参数对精密可控震源压重振动情况影响规律时，考虑大地对震源振动性能的影响，寻找到符合实际的最优参数，从而减小压重振动对震源输出能量的影响，实现提升精密可控震源下传信号能量和信噪比的目的，对于该震源的推广应用具有十分重要的意义。

1 精密可控震源大地频变动力学系统

为研究隔振系统参数对压重振动幅值的影响规律，建立精密可控震源系统动力学模型并进行相关分析计算。

1.1 精密可控震源工作原理

CASS-10国产移动式精密可控震源实物见图 1，它主要由压重、立柱、电机、偏心质量块、激振箱体以及平板组成，2个偏心质量块分别由2个电机直接带动作同步反向旋转，从而产生周期变化的激振力^[14]。

压重提供的重力作用于平板，保证精密可控震源在受到由偏心质量块旋转产生的周期变化的激振力作用时仍与地面保持良好的耦合，避免输出信号发生畸变，且压重质量可根据震源出力要求在一定范围内调节。隔振系统位于立柱和平板之间，一般采用矩形橡胶隔振器。

1.2 基于半无限空间的精密可控震源大地频变动力学模型

根据精密可控震源激振作业的工作原理，可将它转化为两自由度的动力学模型，如图 2所示，图中Q(t)为激振力，m₁为压重的质量，m₂为平板以及激振箱体质量，k₀，c₀是压重与振动平板之间隔振系统的刚度及阻尼，k_d，c_d为耦合地面的刚度及阻尼。

振动器产生的激振力作用于大地表面，而现有研究常将基础当作刚体^{[9, 11, 13]}，作这样假设的计算结果往往与实际有较大的偏差。如图 3所示，实际激振时震源与大地相互作用，大地会影响可控震源的动力学性能，因此基于半无限空间理论开展精密可控震源-大地系统动力学分析。

为简化数学模型，对平板和大地相互作用力进行合理假设：平板变形很小，振动过程中平板与大地不发生脱耦。此时振动器平板对大地的作用力p₀e^iωt近似为均布于半径为r₀的圆面积内的周期力，如图 4所示。

根据半无限空间理论，大地表面受周期扰力p₀e^iωt作用下，Sung以基底中心r=0处的位移w(r, 0, t)近似作为基底位移, 得到了如下结论^[14]：

$ w\left( t \right) = - \frac{{{p_0}{{\rm{e}}^{{\rm{i}}\omega t}}}}{{G{r_0}}}\left( {{f_1} + {\rm{i}}{f_2}} \right) $

(1)

式中：

$ {f_1} = - 0.238\;733 + 0.059\;683a_0^2 - 0.004\;163a_0^4 + \cdots $

$ {f_2} = 0.148\;594{a_0} - 0.017\;757a_0^3 + 0.000\;808a_0^5 - \cdots $

$ {a_0} = \omega {r_0}\sqrt {\rho /G} $

$ G = E/2\left( {1 + \mu } \right) $

其中:ρ, E, μ分别为大地质量密度、弹性模量和泊松比；f₁和f₂为无量纲频率a₀的函数，且随着力的分布不同而变化。

将图 3所示的振动模型中振动装置部分单独分离出来进行分析，如图 5所示，其中R₀e^iωt为大地作用于精密可控震源的反力，Q(t)为激振力。

据公式(1)可知，大地波动反力为：

$ R = {R_0}{{\rm{e}}^{{\rm{i}}\omega t}} = {p_0}{{\rm{e}}^{{\rm{i}}\omega t}} = - \frac{1}{{\left( {{f_1} + {\rm{i}}{f_2}} \right)}}G{r_0}w\left( t \right) $

(2)

对振动装置部分建立动平衡方程：

$ m\ddot \omega \left( t \right) + {R_0}{{\rm{e}}^{{\rm{i}}\omega t}} = Q\left( t \right) $

(3)

式中:m为平板质量，ω(t)为块式平板的位移响应。

将方程(2)代入式(3)即有：

$ m\ddot \omega \left( t \right) + \left[ { - \frac{1}{{\left( {{f_1} + {\rm{i}}{f_2}} \right)}}G{r_0}w\left( t \right)} \right] = Q\left( t \right) $

(4)

对于谐和振动ω(t)，由于iωw(t)=$\dot \omega $(t)且$\omega = \frac{{{a_0}}}{{{r_0}}}\sqrt {G/\rho } $，则式(4)可改写为：

$ \begin{array}{l} m\ddot \omega \left( t \right) + \left[ {\frac{{{f_2}}}{{\left( {f_1^2 + f_2^2} \right){a_0}}}\sqrt {G\rho } r_0^2} \right]\dot \omega \left( t \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\left( { - \frac{{{f_1}}}{{f_1^2 + f_2^2}}G{r_0}} \right)\omega \left( t \right) = Q\left( t \right) \end{array} $

(5)

令${F_1} =-\frac{{{f_1}}}{{f_1^2 + f_2^2}}$, ${F_2} = \frac{{{f_2}}}{{(f_1^2 + f_2^2){a_0}}}$，则式(5)可改为：

$ m\ddot \omega \left( t \right) + \left( {\sqrt {G\rho } r_0^2{F_2}} \right)\dot \omega \left( t \right) + \left( {G{r_0}{F_1}} \right)\omega \left( t \right) = Q\left( t \right) $

(6)

式(6)即是按半无限空间理论模型建立的平板竖直振动形式解。可知半无限空间理论模型中对应耦合地面的刚度系数和阻尼系数表达式分别为：${k_{\rm{d}}} = \sqrt {G\rho } r_0^2{F_2}$，${c_{\rm{d}}} = G{r_0}{F_1}$，且都随频率ω变化。

将k_d，c_d的表达式代入图 2所示的动力学模型中，根据达朗贝尔原理，得到精密可控震源-大地系统频变动力学方程：

$ \mathit{\boldsymbol{M}}\left\{ {\mathit{\boldsymbol{\ddot x}}} \right\} + \mathit{\boldsymbol{C}}\left\{ {\mathit{\boldsymbol{\dot x}}} \right\} + \mathit{\boldsymbol{K}}\left\{ \mathit{\boldsymbol{x}} \right\} = \mathit{\boldsymbol{Q}} $

(7)

其中：

$ \mathit{\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} {{m_1}}&0\\ 0&{{m_2}} \end{array}} \right] $

$ \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {{c_0}}&{ - {c_0}}\\ { - {c_0}}&{{c_{\rm{d}}} + {c_0}} \end{array}} \right] $

$ \mathit{\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {{k_0}}&{ - {k_0}}\\ { - {k_0}}&{{k_{\rm{d}}} + {k_0}} \end{array}} \right] $

系统外力等于震源激振力，根据精密可控震源实际工作的扫描模式，ω(t)线性扫描信号如图 6所示^[15]。

则1个周期内激振频率表达式为:

$ \omega \left( t \right) = \left\{ \begin{array}{l} {\omega _1} + \frac{{2t\left( {{\omega _2} - {\omega _1}} \right)}}{T},\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le t \le \frac{T}{2}\\ {\omega _2} - \frac{{\left( {2t - T} \right)\left( {{\omega _2} - {\omega _1}} \right)}}{T},\;\;\;\;\;\;\frac{T}{2} < t \le T \end{array} \right. $

式中：ω₁，ω₂分别为扫描频宽的下限值和上限值，T为扫描周期。

据精密可控震源出力原理^[15]，激振力表示为：

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} 0\\ {2mr{\omega ^2}\left( t \right)\cos \left( {\omega \left( t \right) \times t} \right)} \end{array}} \right] $

1.3 压重振动剧烈程度评价指标的建立

虽然K，C均为时间t的函数，式(7)属于参变系统方程^[16]，但在震源扫描模式(图 6)下激振频率变化率很小，仅为0.02 Hz/s，而由式(6)可知，k_d，c_d均为ω的函数，为简化计算，在每个时刻将该方程近似看作不变系统方程, 仍按照线性振动系统进行求解。设方程特解(即系统响应)为：

$ \mathit{\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{B_1}{{\rm{e}}^{{\rm{i}}\omega t}}}\\ {{B_2}{{\rm{e}}^{{\rm{i}}\omega t}}} \end{array}} \right] $

其中B₁，B₂为平板与压重的振动幅值，则易得：

$ \mathit{\boldsymbol{\dot x}} = {\rm{i}}\omega \mathit{\boldsymbol{x}},\mathit{\boldsymbol{\ddot x}} = - {\omega ^2}\mathit{\boldsymbol{x}} $

同理可得：

$ \mathit{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} 0\\ F \end{array}} \right]{{\rm{e}}^{{\rm{i}}\omega t}} $

式中F=2mrω²(t)=Aω²(t), 为激振力幅值范围。

将各表达式代入式(7)并化简可得：

$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {{B_1}}\\ {{B_2}} \end{array}} \right] = {\mathit{\boldsymbol{Z}}^{ - 1}}\left[ {\begin{array}{*{20}{c}} 0\\ F \end{array}} \right] = \frac{F}{{\left| \mathit{\boldsymbol{Z}} \right|}}\left[ \begin{array}{l} {k_0} + {\rm{i}}\omega {c_0}\\ {k_0} + {\rm{i}}\omega {c_0} - {m_1}{\omega ^2} \end{array} \right] $

其中：

$ \mathit{\boldsymbol{Z}} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} {k_0} + {\rm{i}}\omega {c_0} - {m_1}{\omega ^2}\\ - {k_0} - {\rm{i}}\omega {c_0} \end{array}&\begin{array}{l} - {k_0} - {\rm{i}}\omega {c_0}\\ {k_{\rm{d}}} + {k_0} + {\rm{i}}\omega {c_{\rm{d}}} + {\rm{i}}\omega {c_{\rm{0}}} - {m_2}{\omega ^2} \end{array} \end{array}} \right] $

由于B₁，B₂均含复数，对它取模即得到精密可控震源压重及平板的幅值表达式。

为了更好地研究精密可控震源压重振动幅值，引入平板受静力下的位移B₀=F/k_d，令：

$ \begin{array}{l} {\beta _1} = {B_1}/{B_0} = {\left[ {{\alpha ^4} + {{\left( {2{\xi _0}\alpha \lambda } \right)}^2}} \right]^{\frac{1}{2}}} \times \left\{ {\left[ {\left( {{\alpha ^2} - {\lambda ^2}} \right)} \right.} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left. {\left( {1 - {\lambda ^2}} \right) - {\mu _{\rm{r}}}{\lambda ^2}{\alpha ^2} - 4{\xi _0}{\xi _{\rm{d}}}\alpha {\lambda ^2}} \right]^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left[ {2{\xi _{\rm{d}}}\lambda \left( {{\alpha ^2} - {\lambda ^2}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left. {{{\left. {2{\xi _0}\alpha \lambda \times \left( {1 - {\mu _{\rm{r}}}{\lambda ^2} - {\lambda ^2}} \right)} \right]}^2}} \right\}^{ - \frac{1}{2}}} \end{array} $

式中：${\xi _0} = \frac{{{c_0}}}{{2{m_1}{\omega _1}}}$，${\xi _{\rm{d}}} = \frac{{{c_{\rm{d}}}}}{{2{m_2}{\omega _2}}}$，$\alpha = \frac{{{\omega _1}}}{{{\omega _2}}}$，$\lambda = \frac{\omega }{{{\omega _2}}}$，${\mu _{\rm{r}}} = \frac{{{m_1}}}{{{m_2}}}$，${\omega _1} = \sqrt {\frac{{{k_0}}}{{{m_1}}}} $，${\omega _2} = \sqrt {\frac{{{k_{\rm{d}}}}}{{{m_2}}}} $。

β₁为压重振幅放大系数，其物理意义为任意时刻压重振幅与平板静位移的比值，该值越大则表示系统对压重振幅放大作用越明显，本文将它作为压重振动剧烈程度的评价指标。

2 隔振参数对评价指标影响规律分析及优化计算

为减少压重对精密可控震源输出信号的影响，需要使振幅放大系数β₁尽可能减小。为此分析隔振参数与β₁之间的关系，以寻找最优的参数组合。精密可控震源与大地各参数如表 1、表 2所示。

将各参数代入β₁的表达式，通过MATLAB软件计算震源扫描频率范围内压重振幅放大系数β₁随各参数的变化规律。

2.1 隔振系统刚度k₀对β₁影响规律分析

压重质量m₁需根据出力以2 000 kg为递增单元在2 000~8 000 kg范围内进行调节。分别计算4 000，6 000 kg压重质量下β₁随k₀的变化情况，结果如图 7所示，其中刚度变化范围参考橡胶隔振器的特性，取k₀=0.1~50 kN·mm^-1[17]。

由图 7(a)可知：1)不同激振频率下压重振幅放大系数极大值所对应的刚度数值并不相同，但均处于5~15 kN·mm^-1这样一个较小的区间内，且频率越高，β₁的极大值时所对应的刚度越大；2)不同激振频率下压重振幅放大系数随隔振系统刚度变化趋势均为先迅速上升后缓慢下降，极小值为0.05，极大值为8.42；3)随着激振频率增大，压重振幅放大系数变化速度变缓。

对比图 7(a)和(b)可知，当压重质量m₁变化时，β₁随k₀的变化规律同样发生改变：1)出现压重振幅放大系数极大值的刚度区间扩大为5~22 kN·mm^-1；2)β₁随刚度变化上升和下降的趋势均变缓。

2.2 隔振系统阻尼比ξ₀对β₁影响规律分析

为了同时表征刚度的影响，分别计算10，40 kN·mm^-1刚度下压重振幅放大系数β₁随阻尼比ξ₀的变化情况，压重质量m₁=6 000 kg，结果如图 8所示，其中阻尼比ξ₀=0.025~0.15^[17]。

由图 8(a)可知：1)随着阻尼比的增大，β₁极大值逐渐减小，且减小趋势逐渐放缓; 2)β₁随着激振频率的增加先上升后下降，在6 Hz处出现极大值，为20.4，在10 Hz处出现极小值，为0.4。

对比图 8(a)和图 8(b)可知，当系统刚度k₀变化时，β₁随ξ₀的变化规律同样发生改变：1)β₁极大值出现在10 Hz处，极大值减小为5.96；2)影响β₁的因素主要为激振频率，而阻尼比对它的影响不明显。

2.3 不同压重质量下隔振系统最优参数组合

由上述分析可知，精密可控震源压重和平板之间的隔振系统的刚度、阻尼比对振幅影响很大，且随着压重质量、激振频率变化时所呈现的动力学规律也有着很大区别。

为得到精密可控震源隔振系统的最优参数组合，基于坐标轮换法设计了优化算法^[18]，其中设计变量为隔振系统刚度k₀和阻尼比ξ₀，以整个频宽范围内平均振幅放大系数β₁最小为优化目标，兼顾系统动力学性能的稳定性和隔振系统的经济性，设立约束条件为：

$ \left\{ \begin{array}{l} 0.1 \le {k_0} \le 50\;{\rm{kN}} \cdot {\rm{m}}{{\rm{m}}^{ - 1}}\\ 0.025 \le {\xi _0} \le 0.15\\ \left| {\frac{{{{\bar \beta }_1}\left( {{{k'}_0} + 1\% \Delta {k_0},{{\xi '}_0}} \right) - {{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}}{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}}} \right| \le 5\% \\ \left| {\frac{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0} + 1\% \Delta {\xi _0}} \right) - {{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}}{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}}} \right| \le 5\% \\ \frac{{{{\bar \beta }_1}\left( {{{k'}_0} + 1\% \Delta {k_0},{{\xi '}_0}} \right)}}{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}} \le 99\% \\ \frac{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0} + 1\% \Delta {{\xi '}_0}} \right)}}{{{{\bar \beta }_1}\left( {{{k'}_0},{{\xi '}_0}} \right)}} \le 99\% \end{array} \right. $

利用MATLAB进行优化计算，其计算逻辑如图 9所示，优化后结果见表 3。

为检验优化效果，以m₁=6 000 kg为例，求解方程(7)，绘制震源扫描信号上半周期(下半周期规律相似)压重振幅的包络线，如图 10所示。

从图中可知，改进前隔振系统在振动过程中压重振幅发生较大的变化，有明显的共振现象发生；改进后压重在震源扫描信号上半周期内振幅变化平缓，且压重振幅极大值较改进前减小了75%。

3 结论

1) 精密可控震源振动动力学性能不仅与自身机械系统参数有关，震源激振频率、近地表地层的弹性模量以及泊松比等物理参数也对它有较大影响。

2) 计算表明不合理的隔振系统参数组合会导致压重发生剧烈振动，整个扫描周期内的压重质量、隔振系统刚度以及阻尼比均对压重振动幅值有较大影响，通过合理设计能有效避免共振现象发生。

3) 通过优化计算得到不同压重质量下，精密可控震源在4~10 Hz扫描频宽范围内的最优刚度及阻尼比参考值，可减小激振作业时压重的振动幅值，以达到削减输出信号干扰、提升与地层之间的能量传递率的目的。

4) 在压重质量为6 000 kg时，经过对隔振参数的优化计算，整个激振周期内优化后的压重振动幅值变化较优化前更加平缓，且振幅极大值下降了75%。