一种海上地震勘探气枪震源子波模拟方法及系统

摘要

本发明提供一种海上地震勘探气枪震源子波模拟方法及系统：测量海水及气枪装置的基本参数；计算：枪体节流常数，气枪枪室的初始温度，枪室中气体的初始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质量，初始静水压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径；确定气泡振动过程中气泡内温度随时间的变化函数，气泡内气体的摩尔质量随时间的变化函数，气泡内的气体体积随时间的变化函数；模拟随时间变化的海水粘度；考虑海水粘度和虚反射压力场，模拟气泡产生的压力波；由气泡产生的压力波，模拟其虚反射压力波；由气泡产生的压力波及虚反射压力波叠加合成气枪子波。本发明增强了模拟子波与实测子波的一致性。

说明书

技术领域

本发明涉及油气地球物理勘探中的地震资料采集领域，尤其涉及一种海上地震勘探气 枪震源子波模拟方法及系统。

背景技术

海上地震勘探是一种利用人工地震技术探测地下构造的勘探方法。它按照一定的方式 在海面附近激发地震波，产生称之为地震子波的振动信号，地震子波由震源开始向下传播， 遇到地质界面之后，在界面处发生透射和反射，透射的地震子波继续向下传播，而反射子 波在界面位置向上传播。来自不同深度界面的地震子波以不同的时间到达海水表面，通过 布设在海面的一种称之为检波器的接收装置，接收来自不同深度地质界面的反射波，接收 到的数字信号称为地震记录。

海上地震采集技术来源于陆上地震勘探，因此，早期的海上勘探也使用陆地勘探的炸 药震源。但很快暴露出其致命弱点：其一，自动化程度差，人工操作危险性大；其二，有 悖于环保理念，对海洋及其海洋生物危害很大。有鉴于此，空气枪、蒸汽枪、水枪等非炸 药震源应运而生。其中，空气枪以其性能稳定、自动化程度高、成本较低等诸多优点逐渐 占据主导地位。

气枪震源所激发的地震子波是海上地震勘探的基本地震信号，其形态特征和频率特征 直接关系到海上地震勘探的最终质量，因此，气枪震源子波数值模拟和特征分析对海上地 震勘探具有十分重要的作用。

Keller和Kolodner（1956）认为气枪产生的高压气体是一个自由振荡的球形气泡，并 提出了“自由气泡振荡原理”。基于这一理论，Ziolkowski（1970）提出了用于模拟气枪 子波的基本数学模型，但是与实测子波相比，Ziolkowski模型所模拟的子波初始脉冲值太 大，而且气泡振荡衰减太慢。Schulze-Gattermann（1972）提出了假设枪体为刚性球体的新 的气泡震荡模型；Ziolkowski（1984）对枪体节流作用和气泡壁热传导作用做出了详细分 析，总结出气体在枪室中并非瞬间释放，而且气体的释放过程会影响气枪子波的形状； Langhammer（1993）通过实验性研究量化了海水粘度对气枪子波的影响。李国发等（2010） 综合考虑枪体节流作用、气泡壁的热传导作用、气泡上浮以及流体粘度等因素对气枪子波 的影响，完成了考虑多种影响因素的单枪子波模拟。

然而，以上模型存在两个基本缺陷，首先，上述模型假设流体粘度为一常数，忽略了 气泡在振动过程中温度和压力的变化对流体粘度的影响。再者，上述模型虽然也考虑了虚 反射对远场子波的影响，但只是简单地将虚反射叠加在气泡振动之上，忽略了虚反射产生 的压力场对气泡振动的再次影响。

发明内容

本发明实施例提供一种海上地震勘探气枪震源子波模拟方法及系统，以建立一种更加 完善的气枪子波模拟系统，使得其所模拟的子波与野外实测子波具有更高的一致性。

一方面，本发明实施例提供了一种海上地震勘探气枪震源子波模拟方法，所述海上地 震勘探气枪震源子波模拟方法包括：

（1）测量并获取海水及其气枪装置的基本参数；

（2）利用所述基本参数计算：枪体节流常数，气枪枪室的初始温度，枪室中气体的 初始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质量，初始静 水压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径；

（3）按照Ziolkowski气泡振动理论，确定气泡振动过程中气泡内温度随时间的变化 函数T(iΔt),i＝0,1,…n，气泡内气体的摩尔质量随时间的变化函数m(iΔt),i＝0,1,…n，气 泡内的气体体积随时间的变化函数V(iΔt),i＝0,1,…n，其中，i时间是采样序号，Δt是采 样间隔，单位是ms,n是时间采样个数；

（4）由气泡内的温度T(iΔt),i＝0,1,…n，模拟随时间变化的海水粘度 μ(iΔt),i＝0,1…n，单位为Kg/(m·s)；

μ(iΔt)＝100+10¹²×T^-6.08(iΔt)；

（5）考虑海水粘度和虚反射压力场对气泡振动的影响，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n；

（6）由气泡产生的压力波w(iΔt),i＝0,1,…n，模拟其虚反射压力波 i＝0,1,…n，单位为Pa，

$w_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0& \mathrm{i\Delta t}<\frac{2h}{c}\\ R\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{c}\end{array}\right)$

其中，R为海水的反射系数，R＝-1，无量纲，h为气枪沉放深度，单位为m，c为 海水中的声波速度，单位为m/ms，c＝1.483m/ms；

（7）由气泡产生的压力波w(iΔt),i＝0,1,…n及虚反射压力波i＝0,1,…n叠 加合成气枪子波W(iΔt),i＝0,1,…n，单位Bar·m。

可选的，在本发明一实施例中，所述（1）中，测量并获取海水及其气枪装置的基本 参数，包括：

（1.1）测量海水温度T_w，单位为K，和海水密度ρ，单位为Kg/m³;

（1.2）测量气枪装置的下列参数：

a.枪室的工作压强p_g，单位为psi；

b.气枪枪室的容积V_g，单位为in³；

c.枪体的释放效率η，无量纲；

d.与枪体容量无关的常数τ₀和节流指数β，无量纲；

e.气枪沉放深度h，单位为m。

可选的，在本发明一实施例中，所述（2）中，利用所述基本参数计算：枪体节流常 数，气枪枪室的初始温度，枪室中气体的初始摩尔质量，气泡内气体的初始摩尔质量，气 泡平衡时气泡内气体的摩尔质量，初始静水压强，气泡的初始半径，气泡中气体的初始体 积，气泡平衡半径，包括：

（2.1）由下式计算枪体节流常数τ，无量纲；

τ＝τ₀(V_g)^β；

（2.2）由海水温度T_w和气枪枪室的初始压强p_g，计算气枪枪室的初始温度T_g，单位 为K，

$T_g=T_w(1+\frac{p_g}{p_c}),$

其中常量p_c＝139Mpa；

（2.3）由枪室初始温度T_g和枪体的释放效率η，计算枪室中初始气体的摩尔质量m_g和 平衡时气泡内的摩尔质量m_g′，单位均为mol，

$m_g=\frac{p_g{V}_g}{T_g{R}_G},$

其中R_G为普氏常量，其值为8.2；

m′_g＝m_gη；

（2.4）由海水密度ρ和气枪沉放深度h，计算初始静水压强p_∞，单位为Pa，

p_∞＝p₀+ρgh，

其中大气压强p₀＝101300Pa，重力加速度g＝9.8m/s²；

（2.5）由初始静水压强p_∞，计算气泡内气体的初始摩尔质量m(0)，单位为mol；

$m\left(0\right)=\frac{p_{\infty }{V}_g}{T_w{R}_G};$

（2.6）由平衡时气泡内气体的摩尔质量m_g′、海水温度T_w和初始静水压强p_∞，计算气 泡的平衡半径r₀以及枪体半径r_g，单位均为m；

$r_0={\left(\frac{3{m_g}^{\prime }{R}_G{T}_w}{4\pi p_{\infty }}\right)}^{\frac{1}{3}},$

r_g＝r₀×0.23；

（2.7）由枪体容积V_g和枪体半径r_g，计算气泡的初始半径R(0)，单位为m；

$R\left(0\right)={\left[\frac{3(V_g+\frac{4}{3}\pi {r_g}^3)}{4\pi }\right]}^{\frac{1}{3}};$

（2.8）由气泡初始半径R(0)和枪体半径r_g计算气泡中气体的初始体积V(0)，单位为 m³；

$V\left(0\right)=\frac{4}{3}\mathrm{\pi R}{\left(0\right)}^3-\frac{4}{3}\pi {r_g}^3.$

可选的，在本发明一实施例中，所述（5）中，考虑海水粘度和虚反射压力场对气泡 振动的影响，模拟气泡产生的压力波w(iΔt),i＝0,1,…n，具体步骤包括：

（5.1）由初始时刻i＝0开始，以递归循环方式模拟气泡产生的压力波 w(iΔt),i＝0,1,…n，其中，R(iΔt)是随时间变化的气泡半径，初值为 单位为m，是气泡半径随时间的变化率，即气泡壁的速 度，初值为单位为m/s；

（5.2）由第i时刻气泡内气体的摩尔质量m(iΔt)、气泡内的气体体积V(iΔt)、气泡 内的温度T(iΔt)和海水粘度μ(iΔt)，模拟第i时刻的气泡壁压力p(iΔt)，单位为Pa，

$p\left(\mathrm{i\Delta t}\right)=\frac{m\left(\mathrm{i\Delta t}\right)R_{G}T\left(\mathrm{i\Delta t}\right)}{V\left(\mathrm{i\Delta t}\right)}+P_V-\frac{2\sigma }{R\left(\mathrm{i\Delta t}\right)}-\frac{4\mu \left(\mathrm{i\Delta t}\right)}{R\left(\mathrm{i\Delta t}\right)}\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

其中，水蒸气压力P_V＝2337Pa，气泡表面张力σ＝0.07061N/m；

（5.3）考虑虚反射压力场对气泡振动的影响，由第i时刻的气泡半径R(iΔt)和气泡壁 的速度模拟第i时刻气泡产生的压力波w(iΔt)和虚反射压力场单位 均为Pa，具体过程为：

首先由气泡半径R(iΔt)和气泡壁的速度得到中间函数f(iΔt)；

$f\left(\mathrm{i\Delta t}\right)=R^2\left(\mathrm{i\Delta t}\right)\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

再利用中间函数f(iΔt)及其对时间导数f′(iΔt)，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n和虚反射压力场

$w\left(\mathrm{i\Delta t}\right)=\rho [-(-f^{\prime }\left(\mathrm{i\Delta t}\right))-\frac{1}{2}{(\frac{f^{\prime }\left(\mathrm{i\Delta t}\right)}{c}+f\left(\mathrm{i\Delta t}\right))}^2],$

$p_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0,& \mathrm{i\Delta t}<\frac{2h}{v}\\ -\frac{1}{1+2h}\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{v}\end{array}\right),$

其中，h为气枪沉放深度，单位为m，c为海水中的声波速度，单位为m/ms， c＝1.483m/ms；

（5.4）由以上结果计算i+1时刻的气泡壁加速度i＝0,1,…n，单位为m/s²，

$\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)=\frac{-\frac{1}{\rho }[p_{\infty }+p_g(\mathrm{i\Delta t}-\frac{2h}{c})-p\left(\mathrm{i\Delta t}\right)-\frac{R\left(\mathrm{i\Delta t}\right)}{c}\stackrel{·}{p}\left(\mathrm{i\Delta t}\right)]-\frac{3}{2}{\stackrel{·}{R}}^2\left(\mathrm{i\Delta t}\right)[1-\frac{4\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{3c}]}{R\left(\mathrm{i\Delta t}\right)[1-\frac{2\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{c}]},$

其中是气泡壁压力p(iΔt)的变化率；

（5.5）由i+1时刻的气泡壁的加速度计算该时刻的气泡壁速度 和气泡半径R((i+1)Δt)；

$\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)=\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)+\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t},$

$R\left((i+1)\mathrm{\Delta t}\right)=R\left(\mathrm{i\Delta t}\right)+\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t};$

（5.6）返回步骤（5.2），计算i+1时刻的气泡产生的压力波w[(i+1)Δt]，由此迭代 方式得到所有时刻的压力波w(iΔt),i＝0,1,…n。

另一方面，本发明实施例提供了一种海上地震勘探气枪震源子波模拟系统，所述海上 地震勘探气枪震源子波模拟系统包括：

测量单元，用于测量并获取海水及其气枪装置的基本参数；

计算单元，用于利用所述基本参数计算：枪体节流常数，气枪枪室的初始温度，枪室 中气体的初始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质量， 初始静水压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径；

函数单元，用于按照Ziolkowski气泡振动理论，确定气泡振动过程中气泡内温度随时 间的变化函数T(iΔt),i＝0,1,…n，气泡内气体的摩尔质量随时间的变化函数 m(iΔt),i＝0,1,…n，气泡内的气体体积随时间的变化函数V(iΔt),i＝0,1,…n，其中，i时间 是采样序号，Δt是采样间隔，单位是ms,n是时间采样个数；

模拟单元，用于由气泡内的温度T(iΔt),i＝0,1,…n，模拟随时间变化的海水粘度 μ(iΔt),i＝0,1…n，单位为Kg/(m·s)；

μ(iΔt)＝100+10¹²×T^-6.08(iΔt)；

用于考虑海水粘度和虚反射压力场对气泡振动的影响，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n；

并用于由气泡产生的压力波w(iΔt),i＝0,1,…n，模拟其虚反射压力波 i＝0,1,…n，单位为Pa，

$w_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0& \mathrm{i\Delta t}<\frac{2h}{c}\\ R\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{c}\end{array}\right)$

其中，R为海水的反射系数，R＝-1，无量纲，h为气枪沉放深度，单位为m，c为 海水中的声波速度，单位为m/ms，c＝1.483m/ms

合成单元，用于由气泡产生的压力波w(iΔt),i＝0,1,…n及虚反射压力波 i＝0,1,…n叠加合成气枪子波W(iΔt),i＝0,1,…n，单位Bar·m。

可选的，在本发明一实施例中，所述测量单元，进一步用于测量并获取海水及其气枪 装置的基本参数，包括：

（1.1）测量海水温度T_w，单位为K，和海水密度ρ，单位为Kg/m³;

（1.2）测量气枪装置的下列参数：

a.枪室的工作压强p_g，单位为psi；

b.气枪枪室的容积V_g，单位为in³；

c.枪体的释放效率η，无量纲；

d.与枪体容量无关的常数τ₀和节流指数β，无量纲；

e.气枪沉放深度h，单位为m。

可选的，在本发明一实施例中，所述计算单元，进一步用于利用所述基本参数计算： 枪体节流常数，气枪枪室的初始温度，枪室中气体的初始摩尔质量，气泡内气体的初始摩 尔质量，气泡平衡时气泡内气体的摩尔质量，初始静水压强，气泡的初始半径，气泡中气 体的初始体积，气泡平衡半径，包括：

（2.1）由下式计算枪体节流常数τ，无量纲；

τ＝τ₀(V_g)^β；

（2.2）由海水温度T_w和气枪枪室的初始压强p_g，计算气枪枪室的初始温度T_g，单位 为K，

$T_g=T_w(1+\frac{p_g}{p_c}),$

其中常量p_c＝139Mpa；

（2.3）由枪室初始温度T_g和枪体的释放效率η，计算枪室中初始气体的摩尔质量m_g和 平衡时气泡内的摩尔质量m_g′，单位均为mol，

$m_g=\frac{p_g{V}_g}{T_g{R}_G},$

其中R_G为普氏常量，其值为8.2；

m_g′＝m_gη；

（2.4）由海水密度ρ和气枪沉放深度h，计算初始静水压强p_∞，单位为Pa，

p_∞＝p₀+ρgh，

其中大气压强p₀＝101300Pa，重力加速度g＝9.8m/s²；

（2.5）由初始静水压强p_∞，计算气泡内气体的初始摩尔质量m(0)，单位为mol；

$m\left(0\right)=\frac{p_{\infty }{V}_g}{T_w{R}_G};$

（2.6）由平衡时气泡内气体的摩尔质量m_g′、海水温度T_w和初始静水压强p_∞，计算气 泡的平衡半径r₀以及枪体半径r_g，单位均为m；

$r_0={\left(\frac{3{m_g}^{\prime }{R}_G{T}_w}{4\pi p_{\infty }}\right)}^{\frac{1}{3}},$

r_g＝r₀×0.23；

（2.7）由枪体容积V_g和枪体半径r_g，计算气泡的初始半径R(0)，单位为m；

$R\left(0\right)={\left[\frac{3(V_g+\frac{4}{3}\pi {r_g}^3)}{4\pi }\right]}^{\frac{1}{3}};$

（2.8）由气泡初始半径R(0)和枪体半径r_g计算气泡中气体的初始体积V(0)，单位为 m³；

$V\left(0\right)=\frac{4}{3}\mathrm{\pi R}{\left(0\right)}^3-\frac{4}{3}\pi {r_g}^3.$

可选的，在本发明一实施例中，所述模拟单元，进一步用于考虑海水粘度和虚反射压 力场对气泡振动的影响，模拟气泡产生的压力波w(iΔt),i＝0,1,…n，具体包括：

（5.1）由初始时刻i＝0开始，以递归循环方式模拟气泡产生的压力波 w(iΔt),i＝0,1,…n，其中，R(iΔt)是随时间变化的气泡半径，初值为 单位为m，是气泡半径随时间的变化率，即气泡壁的速 度，初值为单位为m/s；

（5.2）由第i时刻气泡内气体的摩尔质量m(iΔt)、气泡内的气体体积V(iΔt)、气泡 内的温度T(iΔt)和海水粘度μ(iΔt)，模拟第i时刻的气泡壁压力p(iΔt)，单位为Pa，

$p\left(\mathrm{i\Delta t}\right)=\frac{m\left(\mathrm{i\Delta t}\right)R_{G}T\left(\mathrm{i\Delta t}\right)}{V\left(\mathrm{i\Delta t}\right)}+P_V-\frac{2\sigma }{R\left(\mathrm{i\Delta t}\right)}-\frac{4\mu \left(\mathrm{i\Delta t}\right)}{R\left(\mathrm{i\Delta t}\right)}\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

其中，水蒸气压力P_V＝2337Pa，气泡表面张力σ＝0.07061N/m；

（5.3）考虑虚反射压力场对气泡振动的影响，由第i时刻的气泡半径R(iΔt)和气泡壁 的速度模拟第i时刻气泡产生的压力波w(iΔt)和虚反射压力场单位 均为Pa，具体过程为：

首先由气泡半径R(iΔt)和气泡壁的速度得到中间函数f(iΔt)；

$f\left(\mathrm{i\Delta t}\right)=R^2\left(\mathrm{i\Delta t}\right)\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

再利用中间函数f(iΔt)及其对时间的导数f′(iΔt)，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n和虚反射压力场

$w\left(\mathrm{i\Delta t}\right)=\rho [-(-f^{\prime }\left(\mathrm{i\Delta t}\right))-\frac{1}{2}{(\frac{f^{\prime }\left(\mathrm{i\Delta t}\right)}{c}+f\left(\mathrm{i\Delta t}\right))}^2],$

$p_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0,& \mathrm{i\Delta t}<\frac{2h}{v}\\ -\frac{1}{1+2h}\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{v}\end{array}\right),$

其中，h为气枪沉放深度，单位为m，c为海水中的声波速度，单位为m/ms， c＝1.483m/ms；

（5.4）由以上结果计算i+1时刻的气泡壁加速度i＝0,1,…n，单位为m/s²，

$\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)=\frac{-\frac{1}{\rho }[p_{\infty }+p_g(\mathrm{i\Delta t}-\frac{2h}{c})-p\left(\mathrm{i\Delta t}\right)-\frac{R\left(\mathrm{i\Delta t}\right)}{c}\stackrel{·}{p}\left(\mathrm{i\Delta t}\right)]-\frac{3}{2}{\stackrel{·}{R}}^2\left(\mathrm{i\Delta t}\right)[1-\frac{4\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{3c}]}{R\left(\mathrm{i\Delta t}\right)[1-\frac{2\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{c}]},$

其中是气泡壁压力p(iΔt)的变化率；

（5.5）由i+1时刻的气泡壁的加速度计算该时刻的气泡壁速度 和气泡半径R((i+1)Δt)；

$\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)=\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)+\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t},$

$R\left((i+1)\mathrm{\Delta t}\right)=R\left(\mathrm{i\Delta t}\right)+\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t};$

（5.6）返回步骤（5.2），计算i+1时刻的气泡产生的压力波w[(i+1)Δt]，由此迭代 方式得到所有时刻的压力波w(iΔt),i＝0,1,…n

上述技术方案具有如下有益效果：因为采用综合考虑流体粘度随温度变化以及虚反射 对气泡周围压力的影响的技术手段，所以达到了如下的技术效果：建立了一种更加完善的 气枪子波模拟系统，增强模拟子波与实测子波的一致性，为海上地震勘探提供基础数据， 提高海上地震勘探精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技 术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明 的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根 据这些附图获得其他的附图。

图1为本发明实施例一种海上地震勘探气枪震源子波模拟方法流程图；

图2为本发明实施例一种海上地震勘探气枪震源子波模拟系统结构示意图；

图3为本发明应用实例本系统所模拟的气泡内温度随时间的变化；

图4为本发明应用实例本系统所模拟气泡中气体质量随时间的变化；

图5为本发明应用实例本系统所模拟气泡中气体体积随时间的变化；

图6为本发明应用实例本系统所模拟海水粘度随温度的变化；

图7为本发明应用实例本系统所模拟的气泡产生的压力波随时间的变化；

图8为本发明应用实例本系统所模拟的气泡半径随时间的变化；

图9为本发明应用实例本系统所模拟的气泡壁速度随时间的变化；

图10为本发明应用实例本系统所模拟的气泡壁加速度随时间的变化；

图11为本发明应用实例本系统所模拟的气泡壁压力随时间的变化；

图12为本发明应用实例本系统所模拟的气泡虚反射压力波随时间的变化；

图13为本发明应用实例本系统所模拟的气枪子波随时间的变化；

图14为本发明应用实例某海上工区实际激发的气枪子波与按照李国发等所提供的方 法模拟的气枪子波对比示意图；

图15为本发明应用实例某海上工区实际激发的气枪子波与本发明所模拟的气枪子波 对比示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地 描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本 发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实 施例，都属于本发明保护的范围。

本发明实施例提供一种海上地震勘探气枪震源子波模拟方法及系统，以建立一种更加 完善的气枪子波模拟系统，使得其所模拟的子波与野外实测子波具有更高的一致性。

如图1所示，为本发明实施例一种海上地震勘探气枪震源子波模拟方法流程图，所述 海上地震勘探气枪震源子波模拟方法包括：

（1）测量并获取海水及其气枪装置的基本参数；

（2）利用所述基本参数计算：枪体节流常数，气枪枪室的初始温度，枪室中气体的 初始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质量，初始静 水压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径；

（3）按照Ziolkowski气泡振动理论，确定气泡振动过程中气泡内温度随时间的变化 函数T(iΔt),i＝0,1,…n，气泡内气体的摩尔质量随时间的变化函数m(iΔt),i＝0,1,…n，气 泡内的气体体积随时间的变化函数V(iΔt),i＝0,1,…n，其中，i时间是采样序号，Δt是采 样间隔，单位是ms,n是时间采样个数；

（4）由气泡内的温度T(iΔt),i＝0,1,…n，模拟随时间变化的海水粘度 μ(iΔt),i＝0,1…n，单位为Kg/(m·s)；μ(iΔt)＝100+10¹²×T^-6.08(iΔt)；

（5）考虑海水粘度和虚反射压力场对气泡振动的影响，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n；

（6）由气泡产生的压力波w(iΔt),i＝0,1,…n，模拟其虚反射压力波 i＝0,1,…n，单位为Pa，

$w_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0& \mathrm{i\Delta t}<\frac{2h}{c}\\ R\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{c}\end{array}\right)$

其中，R为海水的反射系数，R＝-1，无量纲，h为气枪沉放深度，单位为m，c为 海水中的声波速度，单位为m/ms，c＝1.483m/ms；

（7）由气泡产生的压力波w(iΔt),i＝0,1,…n及虚反射压力波i＝0,1,…n叠 加合成气枪子波W(iΔt),i＝0,1,…n，单位Bar·m。

可选的，在本发明一实施例中，所述（1）中，测量并获取海水及其气枪装置的基本 参数，包括：

（1.1）测量海水温度T_w，单位为K，和海水密度ρ，单位为Kg/m³;

（1.2）测量气枪装置的下列参数：

a.枪室的工作压强p_g，单位为psi；

b.气枪枪室的容积V_g，单位为in³；

c.枪体的释放效率η，无量纲；

d.与枪体容量无关的常数τ₀和节流指数β，无量纲；

e.气枪沉放深度h，单位为m。

可选的，在本发明一实施例中，所述（2）中，利用所述基本参数计算：枪体节流常 数，气枪枪室的初始温度，枪室中气体的初始摩尔质量，气泡内气体的初始摩尔质量，气 泡平衡时气泡内气体的摩尔质量，初始静水压强，气泡的初始半径，气泡中气体的初始体 积，气泡平衡半径，包括：

（2.1）由下式计算枪体节流常数τ，无量纲；

τ＝τ₀(V_g)^β；

（2.2）由海水温度T_w和气枪枪室的初始压强p_g，计算气枪枪室的初始温度T_g，单位 为K，

$T_g=T_w(1+\frac{p_g}{p_c}),$

其中常量p_c＝139Mpa；

（2.3）由枪室初始温度T_g和枪体的释放效率η，计算枪室中初始气体的摩尔质量m_g和 平衡时气泡内的摩尔质量m_g′，单位均为mol，

$m_g=\frac{p_g{V}_g}{T_g{R}_G},$

其中R_G为普氏常量，其值为8.2；

m_g′＝m_gη；

（2.4）由海水密度ρ和气枪沉放深度h，计算初始静水压强p_∞，单位为Pa，

p_∞＝p₀+ρgh，

其中大气压强p₀＝101300Pa，重力加速度g＝9.8m/s²；

（2.5）由初始静水压强p_∞，计算气泡内气体的初始摩尔质量m(0)，单位为mol；

$m\left(0\right)=\frac{p_{\infty }{V}_g}{T_w{R}_G};$

（2.6）由平衡时气泡内气体的摩尔质量m_g′、海水温度T_w和初始静水压强p_∞，计算气 泡的平衡半径r₀以及枪体半径r_g，单位均为m；

$r_0={\left(\frac{3{m_g}^{\prime }{R}_G{T}_w}{4\pi p_{\infty }}\right)}^{\frac{1}{3}},$

r_g＝r₀×0.23；

（2.7）由枪体容积V_g和枪体半径r_g，计算气泡的初始半径R(0)，单位为m；

$R\left(0\right)={\left[\frac{3(V_g+\frac{4}{3}\pi {r_g}^3)}{4\pi }\right]}^{\frac{1}{3}};$

（2.8）由气泡初始半径R(0)和枪体半径r_g计算气泡中气体的初始体积V(0)，单位为 m³；

$V\left(0\right)=\frac{4}{3}\mathrm{\pi R}{\left(0\right)}^3-\frac{4}{3}\pi {r_g}^3.$

可选的，在本发明一实施例中，所述（5）中，考虑海水粘度和虚反射压力场对气泡 振动的影响，模拟气泡产生的压力波w(iΔt),i＝0,1,…n，具体步骤包括：

（5.1）由初始时刻i＝0开始，以递归循环方式模拟气泡产生的压力波 w(iΔt),i＝0,1,…n，其中，R(iΔt)是随时间变化的气泡半径，初值为 单位为m，是气泡半径随时间的变化率，即气泡壁的速 度，初值为单位为m/s；

（5.2）由第i时刻气泡内气体的摩尔质量m(iΔt)、气泡内的气体体积V(iΔt)、气泡 内的温度T(iΔt)和海水粘度μ(iΔt)，模拟第i时刻的气泡壁压力p(iΔt)，单位为Pa，

$p\left(\mathrm{i\Delta t}\right)=\frac{m\left(\mathrm{i\Delta t}\right)R_{G}T\left(\mathrm{i\Delta t}\right)}{V\left(\mathrm{i\Delta t}\right)}+P_V-\frac{2\sigma }{R\left(\mathrm{i\Delta t}\right)}-\frac{4\mu \left(\mathrm{i\Delta t}\right)}{R\left(\mathrm{i\Delta t}\right)}\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

其中，水蒸气压力P_V＝2337Pa，气泡表面张力σ＝0.07061N/m；

（5.3）考虑虚反射压力场对气泡振动的影响，由第i时刻的气泡半径R(iΔt)和气泡壁 的速度模拟第i时刻气泡产生的压力波w(iΔt)和虚反射压力场单位 均为Pa，具体过程为：

首先由气泡半径R(iΔt)和气泡壁的速度得到中间函数f(iΔt)；

$f\left(\mathrm{i\Delta t}\right)=R^2\left(\mathrm{i\Delta t}\right)\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

再利用中间函数f(iΔt)及其对时间的导数f′(iΔt)，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n和虚反射压力场

$w\left(\mathrm{i\Delta t}\right)=\rho [-(-f^{\prime }\left(\mathrm{i\Delta t}\right))-\frac{1}{2}{(\frac{f^{\prime }\left(\mathrm{i\Delta t}\right)}{c}+f\left(\mathrm{i\Delta t}\right))}^2],$

$p_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0,& \mathrm{i\Delta t}<\frac{2h}{v}\\ -\frac{1}{1+2h}\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{v}\end{array}\right),$

其中，h为气枪沉放深度，单位为m，c为海水中的声波速度，单位为m/ms， c＝1.483m/ms；

（5.4）由以上结果计算i+1时刻的气泡壁加速度i＝0,1,…n，单位为m/s²，

$\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)=\frac{-\frac{1}{\rho }[p_{\infty }+p_g(\mathrm{i\Delta t}-\frac{2h}{c})-p\left(\mathrm{i\Delta t}\right)-\frac{R\left(\mathrm{i\Delta t}\right)}{c}\stackrel{·}{p}\left(\mathrm{i\Delta t}\right)]-\frac{3}{2}{\stackrel{·}{R}}^2\left(\mathrm{i\Delta t}\right)[1-\frac{4\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{3c}]}{R\left(\mathrm{i\Delta t}\right)[1-\frac{2\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{c}]},$

其中是气泡壁压力p(iΔt)的变化率；

（5.5）由i+1时刻的气泡壁的加速度计算该时刻的气泡壁速度 和气泡半径R((i+1)Δt)；

$\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)=\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)+\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t},$

$R\left((i+1)\mathrm{\Delta t}\right)=R\left(\mathrm{i\Delta t}\right)+\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t};$

（5.6）返回步骤（5.2），计算i+1时刻的气泡产生的压力波w「(i+1)Δt]由此迭代 方式得到所有时刻的压力波w(iΔt),i＝0,1,…n。

如图2所示，为本发明实施例一种海上地震勘探气枪震源子波模拟系统结构示意图， 所述海上地震勘探气枪震源子波模拟系统包括：

测量单元21，用于测量并获取海水及其气枪装置的基本参数；

计算单元22，用于利用所述基本参数计算：枪体节流常数，气枪枪室的初始温度，枪 室中气体的初始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质 量，初始静水压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径；

函数单元23，用于按照Ziolkowski气泡振动理论，确定气泡振动过程中气泡内温度随 时间的变化函数T(iΔt),i＝0,1,…n，气泡内气体的摩尔质量随时间的变化函数 m(iΔt),i＝0,1,…n，气泡内的气体体积随时间的变化函数V(iΔt),i＝0,1,…n，其中，i时间 是采样序号，Δt是采样间隔，单位是ms,n是时间采样个数；

模拟单元24，用于由气泡内的温度T(iΔt),i＝0,1,…n，模拟随时间变化的海水粘度 μ(iΔt),i＝0,1…n，单位为Kg/(m·s)；

μ(iΔt)＝100+10¹²×T^-6.08(iΔt)；

用于考虑海水粘度和虚反射压力场对气泡振动的影响，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n；

并用于由气泡产生的压力波w(iΔt),i＝0,1,…n，模拟其虚反射压力波 i＝0,1,…n，单位为Pa，

$w_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0& \mathrm{i\Delta t}<\frac{2h}{c}\\ R\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{c}\end{array}\right)$

其中，R为海水的反射系数，R＝-1，无量纲，h为气枪沉放深度，单位为m，c为 海水中的声波速度，单位为m/ms，c＝1.483m/ms；

合成单元25，用于由气泡产生的压力波w(iΔt),i＝0,1,…n及虚反射压力波 i＝0,1,…n叠加合成气枪子波W(iΔt),i＝0,1,…n，单位Bar·m。

可选的，在本发明一实施例中，所述测量单元21，进一步用于测量并获取海水及其气 枪装置的基本参数，包括：

（1.1）测量海水温度T_w，单位为K，和海水密度ρ，单位为Kg/m³;

（1.2）测量气枪装置的下列参数：

a.枪室的工作压强p_g，单位为psi；

b.气枪枪室的容积V_g，单位为in³；

c.枪体的释放效率η，无量纲；

d.与枪体容量无关的常数τ₀和节流指数β，无量纲；

e.气枪沉放深度h，单位为m。

可选的，在本发明一实施例中，所述计算单元22，进一步用于利用所述基本参数计 算：枪体节流常数，气枪枪室的初始温度，枪室中气体的初始摩尔质量，气泡内气体的初 始摩尔质量，气泡平衡时气泡内气体的摩尔质量，初始静水压强，气泡的初始半径，气泡 中气体的初始体积，气泡平衡半径，包括：

（2.1）由下式计算枪体节流常数τ，无量纲；

τ＝τ₀(V_g)^β；

（2.2）由海水温度T_w和气枪枪室的初始压强p_g，计算气枪枪室的初始温度T_g，单位 为K，

其中常量p_c＝139Mpa；

（2.3）由枪室初始温度T_g和枪体的释放效率η，计算枪室中初始气体的摩尔质量m_g和 平衡时气泡内的摩尔质量m_g′，单位均为mol，

$m_g=\frac{p_g{V}_g}{T_g{R}_G},$

其中R_G为普氏常量，其值为8.2；

m_g′＝m_gη；

（2.4）由海水密度ρ和气枪沉放深度h，计算初始静水压强p_∞，单位为Pa，

p_∞＝p₀+ρgh，

其中大气压强p₀＝101300Pa，重力加速度g＝9.8m/s²；

（2.5）由初始静水压强p_∞，计算气泡内气体的初始摩尔质量m(0)，单位为mol；

$m\left(0\right)=\frac{p_{\infty }{V}_g}{T_w{R}_G};$

（2.6）由平衡时气泡内气体的摩尔质量m_g′、海水温度T_w和初始静水压强p_∞，计算气 泡的平衡半径r₀以及枪体半径r_g，单位均为m；

$r_0={\left(\frac{3{m_g}^{\prime }{R}_G{T}_w}{4\pi p_{\infty }}\right)}^{\frac{1}{3}},$

r_g＝r₀×0.23；

（2.7）由枪体容积V_g和枪体半径r_g，计算气泡的初始半径R(0)，单位为m；

$R\left(0\right)={\left[\frac{3(V_g+\frac{4}{3}\pi {r_g}^3)}{4\pi }\right]}^{\frac{1}{3}};$

（2.8）由气泡初始半径R(0)和枪体半径r_g计算气泡中气体的初始体积V(0)，单位为 m³；

$V\left(0\right)=\frac{4}{3}\mathrm{\pi R}{\left(0\right)}^3-\frac{4}{3}\pi {r_g}^3.$

可选的，在本发明一实施例中，所述模拟单元24，进一步用于考虑海水粘度和虚反 射压力场对气泡振动的影响，模拟气泡产生的压力波w(iΔt),i＝0,1,…n，具体包括：

（5.1）由初始时刻i＝0开始，以递归循环方式模拟气泡产生的压力波 w(iΔt),i＝0,1,…n，其中，R(iΔt)是随时间变化的气泡半径，初值为 单位为m，是气泡半径随时间的变化率，即气泡壁的速 度，初值为单位为m/s；

（5.2）由第i时刻气泡内气体的摩尔质量m(iΔt)、气泡内的气体体积V(iΔt)、气泡 内的温度T(iΔt)和海水粘度μ(iΔt)，模拟第i时刻的气泡壁压力p(iΔt)，单位为Pa，

$p\left(\mathrm{i\Delta t}\right)=\frac{m\left(\mathrm{i\Delta t}\right)R_{G}T\left(\mathrm{i\Delta t}\right)}{V\left(\mathrm{i\Delta t}\right)}+P_V-\frac{2\sigma }{R\left(\mathrm{i\Delta t}\right)}-\frac{4\mu \left(\mathrm{i\Delta t}\right)}{R\left(\mathrm{i\Delta t}\right)}\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

其中，水蒸气压力P_V＝2337Pa，气泡表面张力σ＝0.07061N/m；

（5.3）考虑虚反射压力场对气泡振动的影响，由第i时刻的气泡半径R(iΔt)和气泡壁 的速度模拟第i时刻气泡产生的压力波w(iΔt)和虚反射压力场单位 均为Pa，具体过程为：

首先由气泡半径R(iΔt)和气泡壁的速度得到中间函数f(iΔt)；

$f\left(\mathrm{i\Delta t}\right)=R^2\left(\mathrm{i\Delta t}\right)\stackrel{·}{R}\left(\mathrm{i\Delta t}\right),$

再利用中间函数f(iΔt)及其对时间的导数f′(iΔt)，模拟气泡产生的压力波 w(iΔt),i＝0,1,…n和虚反射压力场

$w\left(\mathrm{i\Delta t}\right)=\rho [-(-f^{\prime }\left(\mathrm{i\Delta t}\right))-\frac{1}{2}{(\frac{f^{\prime }\left(\mathrm{i\Delta t}\right)}{c}+f\left(\mathrm{i\Delta t}\right))}^2],$

$p_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0,& \mathrm{i\Delta t}<\frac{2h}{v}\\ -\frac{1}{1+2h}\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{v}\end{array}\right),$

其中，h为气枪沉放深度，单位为m，c为海水中的声波速度，单位为m/ms， c＝1.483m/ms；

（5.4）由以上结果计算i+1时刻的气泡壁加速度i＝0,1,…n，单位为m/s²，

$\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)=\frac{-\frac{1}{\rho }[p_{\infty }+p_g(\mathrm{i\Delta t}-\frac{2h}{c})-p\left(\mathrm{i\Delta t}\right)-\frac{R\left(\mathrm{i\Delta t}\right)}{c}\stackrel{·}{p}\left(\mathrm{i\Delta t}\right)]-\frac{3}{2}{\stackrel{·}{R}}^2\left(\mathrm{i\Delta t}\right)[1-\frac{4\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{3c}]}{R\left(\mathrm{i\Delta t}\right)[1-\frac{2\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{c}]},$

其中是气泡壁压力p(iΔt)的变化率；

（5.5）由i+1时刻的气泡壁的加速度计算该时刻的气泡壁速度 和气泡半径R((i+1)Δt)；

$\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)=\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)+\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t},$

$R\left((i+1)\mathrm{\Delta t}\right)=R\left(\mathrm{i\Delta t}\right)+\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t};$

（5.6）返回步骤（5.2），计算i+1时刻的气泡产生的压力波w[(i+1)Δt]由此迭代 方式得到所有时刻的压力波w(iΔt),i＝0,1,…n。

本发明上述方法或系统实施例上述技术方案具有如下有益效果：因为采用综合考虑流 体粘度随温度变化以及虚反射对气泡周围压力的影响的技术手段，所以达到了如下的技术 效果：建立了一种更加完善的气枪子波模拟系统，增强模拟子波与实测子波的一致性，为 海上地震勘探提供基础数据，提高海上地震勘探精度。

以下举应用实例进行详细说明：

本发明应用实例通过综合考虑流体粘度随温度的变化以及虚反射对气泡周围压力的 影响，建立一种更完善的模拟气枪子波的模型，使其所模拟的子波与野外实测子波具有更 高的一致性。

下面以某气枪震源海上勘探的实际数据为例。如图3所示，为本发明应用实例本系统 所模拟的气泡内温度随时间的变化；如图4所示，为本发明应用实例本系统所模拟气泡中 气体质量随时间的变化；如图5所示，为本发明应用实例本系统所模拟气泡中气体体积随 时间的变化；如图6所示，为本发明应用实例本系统所模拟海水粘度随温度的变化；如图 7所示，为本发明应用实例本系统所模拟的气泡产生的压力波随时间的变化；如图8所示， 为本发明应用实例本系统所模拟的气泡半径随时间的变化；如图9所示，为本发明应用实 例本系统所模拟的气泡壁速度随时间的变化；如图10所示，为本发明应用实例本系统所 模拟的气泡壁加速度随时间的变化；如图11所示，为本发明应用实例本系统所模拟的气 泡壁压力随时间的变化；如图12所示，为本发明应用实例本系统所模拟的气泡虚反射压 力波随时间的变化；如图13所示，为本发明应用实例本系统所模拟的气枪子波随时间的 变化。以下说明具体实施过程：

（1）测量海水及其气枪装置的基本参数：

（1.1）测量海水温度T_w＝300K和海水密度ρ＝1000Kg/m³;

（1.2）测量气枪装置的下列参数：

a.枪室的工作压强p_g＝2000psi；

b.气枪枪室的容积V_g＝40in³；

c.枪体的释放效率η＝0.8317，无量纲；

d.与枪体容量无关的常数τ₀＝0.12和节流指数β＝0.52，无量纲；

e.气枪沉放深度h＝5m；

（2）利用上述测量参数计算枪体节流常数，气枪枪室的初始温度，枪室中气体的初 始摩尔质量，气泡内气体的初始摩尔质量，气泡平衡时气泡内气体的摩尔质量，初始静水 压强，气泡的初始半径，气泡中气体的初始体积，气泡平衡半径。

（2.1）由下式计算枪体节流常数τ，无量纲；

τ＝τ₀(V_g)^β＝0.0027

（2.2）由海水温度T_w和气枪枪室的初始压强p_g，计算气枪枪室的初始温度T_g，单位 为K

$T_g=T_w(1+\frac{p_g}{p_c})=330.2190K$

其中常量p_c＝139Mpa；

（2.3）由枪室初始温度T_g和枪体的释放效率η，计算枪室中初始气体的摩尔质量m_g和 平衡时气泡内的摩尔质量m_g′，单位均为mol

$m_g=\frac{p_g{V}_g}{T_g{R}_G}=3.3406\mathrm{mol}$

其中R_G为普氏常量，其值为8.2；

m_g′＝m_gη＝2.7784mol

（2.4）由海水密度ρ和气枪沉放深度h，计算初始静水压强p_∞，单位为Pa p_∞＝p₀+ρgh＝135900Pa

其中大气压强p₀＝101300Pa，重力加速度g＝9.8m/s²；

（2.5）由初始静水压强p_∞，计算气泡内气体的初始摩尔质量m(0)，单位为mol；

$m\left(0\right)=\frac{p_{\infty }{V}_g}{T_w{R}_G}=0.04\mathrm{mol}$

（2.6）由平衡时气泡内气体的摩尔质量m_g′、海水温度T_w和初始静水压强p_∞，计算气 泡的平衡半径r₀以及枪体半径r_g，单位均为m；

$r_0={\left(\frac{3{m_g}^{\prime }{R}_G{T}_w}{4\pi p_{\infty }}\right)}^{\frac{1}{3}}=0.2214m$

r_g＝r₀×0.23＝0.0509m

（2.7）由枪体容积V_g和枪体半径r_g，计算气泡的初始半径R(0)，单位为m；

$R\left(0\right)={\left[\frac{3(V_g+\frac{4}{3}\pi {r_g}^3)}{4\pi }\right]}^{\frac{1}{3}}=0.0661m$

（2.8）由气泡初始半径R(0)和枪体半径r_g计算气泡中气体的初始体积V(0)，单位为 m³；

$V\left(0\right)=\frac{4}{3}\mathrm{\pi R}{\left(0\right)}^3-\frac{4}{3}\pi {r_g}^3=0.00065548m^3$

（3）按照Ziolkowski气泡振动理论，确定气泡振动过程中如图3所示气泡内温度随 时间的变化函数T(iΔt),i＝0,1,…n，如图4所示的气泡内气体摩尔质量随时间的变化函数 m(iΔt),i＝0,1,…n，如图5所示的气泡内气体体积随时间的变化函数V(iΔt),i＝0,1,…n， 其中，i时间是采样序号，Δt是采样间隔，单位是ms,n是时间采样个数，本例中，采样间 隔Δt＝0.1ms，采样个数n＝3001；

（4）由气泡内的温度T(iΔt),i＝0,1,…n，模拟如图6所示随时间变化的海水粘度 μ(iΔt),i＝0,1…n，单位为Kg/(m·s)；

μ(iΔt)＝100+10¹²×T^-6.08(iΔt)

（5）考虑海水粘度和虚反射压力场对气泡振动的影响，模拟如图7所示气泡产生的 压力波w(iΔt),i＝0,1,…n，具体步骤为：

（5.1）由初始时刻i＝0开始，以递归循环方式模拟气泡产生的压力波 w(iΔt),i＝0,1,…n，其中，R(iΔt)是随时间变化的气泡半径，初值为 是气泡半径随时间的变化率，即气泡壁的速度， 初值为

（5.2）由第i时刻气泡内气体的摩尔质量m(iΔt)、气泡内的气体体积V(iΔt)、气泡内 的温度T(iΔt)和海水粘度μ(iΔt)，模拟第i时刻的气泡壁压力p(iΔt)，单位为Pa

$p\left(\mathrm{i\Delta t}\right)=\frac{m\left(\mathrm{i\Delta t}\right)R_{G}T\left(\mathrm{i\Delta t}\right)}{V\left(\mathrm{i\Delta t}\right)}+P_V-\frac{2\sigma }{R\left(\mathrm{i\Delta t}\right)}-\frac{4\mu \left(\mathrm{i\Delta t}\right)}{R\left(\mathrm{i\Delta t}\right)}\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)$

其中，水蒸气压力P_V＝2337Pa，气泡表面张力σ＝0.07061N/m；

（5.3）考虑虚反射压力场对气泡振动的影响，由第i时刻的气泡半径R(iΔt)和气泡壁 的速度模拟第i时刻气泡产生的压力波w(iΔt)和虚反射压力场单位 均为Pa，具体过程为：

首先由气泡半径R(iΔt)和气泡壁的速度得到中间函数f(iΔt)；

$f\left(\mathrm{i\Delta t}\right)=R^2\left(\mathrm{i\Delta t}\right)\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)$

再利用中间函数f(iΔt)及其对时间的导数f′(iΔt)，模拟气泡产生的压力波w(iΔt) 和虚反射压力场

$w\left(\mathrm{i\Delta t}\right)=\rho [-(-f^{\prime }\left(\mathrm{i\Delta t}\right))-\frac{1}{2}{(\frac{f^{\prime }\left(\mathrm{i\Delta t}\right)}{c}+f\left(\mathrm{i\Delta t}\right))}^2]$

$p_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0,& \mathrm{i\Delta t}<\frac{2h}{v}\\ -\frac{1}{1+2h}\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{v}\end{array}\right)$

其中，h为气枪沉放深度，单位为m，c为海水中的声波速度，单位为m/ms， c＝1.483m/ms；

（5.4）由以上结果计算i+1时刻气泡壁的加速度单位为m/s²

$\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)=\frac{-\frac{1}{\rho }[p_{\infty }+p_g(\mathrm{i\Delta t}-\frac{2h}{c})-p\left(\mathrm{i\Delta t}\right)-\frac{R\left(\mathrm{i\Delta t}\right)}{c}\stackrel{·}{p}\left(\mathrm{i\Delta t}\right)]-\frac{3}{2}{\stackrel{·}{R}}^2\left(\mathrm{i\Delta t}\right)[1-\frac{4\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{3c}]}{R\left(\mathrm{i\Delta t}\right)[1-\frac{2\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)}{c}]}$

其中是气泡壁压力p(iΔt)的变化率；

（5.5）由第i+1时刻的气泡壁加速度计算该时刻的气泡壁速度 和气泡半径R((i+1)Δt)；

$\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)=\stackrel{·}{R}\left(\mathrm{i\Delta t}\right)+\stackrel{··}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t}$

$R\left((i+1)\mathrm{\Delta t}\right)=R\left(\mathrm{i\Delta t}\right)+\stackrel{·}{R}\left((i+1)\mathrm{\Delta t}\right)\mathrm{\Delta t}$

（5.6）返回步骤（5.2），计算i+1时刻气泡产生的压力波w[(i+1)Δt]，由此迭代 方式得到所有时刻的压力波w(iΔt),i＝0,1,…n。在此过程中产生的中间成果有，如图8所 示的气泡半径R(iΔt),i＝0,1,…n随时间的变化情况，如图9所示的气泡半径随时间的变化 率即气泡壁速度i＝0,1,…n随时间的变化情况，如图10所示的气泡壁加速度 i＝0,1,…n随时间的变化情况，如图11所示的气泡壁压力p(iΔt),i＝0,1,…n随 时间的变化情况。

（6）由气泡产生的压力波w(iΔt),i＝0,1,…n，模拟如图12所示虚反射压力波 i＝0,1,…n，单位为Pa

$w_g(\mathrm{i\Delta t}-\frac{2h}{c})=\left(\begin{array}{cc}0& \mathrm{i\Delta t}<\frac{2h}{c}\\ R\times w\left(\mathrm{i\Delta t}\right),& \mathrm{i\Delta t}\ge \frac{2h}{c}\end{array}\right)$

其中，R为海水的反射系数，R＝-1，无量纲，h为气枪沉放深度，单位为m，c为 海水中的声波速度，单位为m/ms，c＝1.483m/ms；

（7）由气泡产生的压力波w(iΔt),i＝0,1,…n及虚反射压力波i＝0,1,…n叠 加合成如图13所示气枪子波W(iΔt),i＝0,1,…n，单位Bar·m。

例如：本实施例为某气枪震源海上地震勘探工区的应用实例，其中气枪容量V_g＝80in³， 陈放深度h＝8m，气枪工作压力p_g＝2000psi；海水中声波速度为c＝1.483m/ms，海水温 度T_w＝20℃。

本实施例气枪子波数据采样间隔为0.1ms,记录长度为300ms。如图14所示，为本发明 应用实例某海上工区实际激发的气枪子波与按照李国发等（2010）所提供的方法（具体方 法参见Li Guofa,Cao Mingqiang,et al.2010,Modeling of the Signature of Air Gun in Marine Seismic Exploration Considering the effects of multiple Practical Physics.Applied Geophysics，7(2):158-165）模拟 的气枪子波对比示意图，其中，实线是该海上工区实际激发的气枪子波，虚线是按照李国 发等（2010）所提供的方法在没有考虑流体粘度随温度的变化及其虚反射压力场对气泡振 动影响情况下模拟的气枪子波。如图15所示，为本发明应用实例某海上工区实际激发的 气枪子波与本发明所模拟的气枪子波对比示意图，其中，实线是该海上工区实际激发的气 枪子波，虚线是本发明的方法在考虑流体粘度随温度的变化及其虚反射压力场对气泡振动 影响情况下模拟的气枪子波。与图14进行对比可以看出，本发明所模拟的子波与实际子 波具有更高的一致性。

本发明综合考虑流体粘度随温度的变化以及虚反射对气泡周围压力的影响，对传统气 枪子波模拟模型进行修正，建立了一种更加完善的气枪子波模拟系统，所模拟的子波与野 外实测子波具有更高的一致性，对海上地震勘探具有重要的指导意义和参考价值。

本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块（illustrative logical block），单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为 清楚展示硬件和软件的可替换性（interchangeability），上述的各种说明性部件（illustrative components），单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软 件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的 应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保 护的范围。

本发明实施例中所描述的各种说明性的逻辑块，或单元都可以通过通用处理器，数字 信号处理器，专用集成电路（ASIC），现场可编程门阵列（FPGA）或其它可编程逻辑装 置，离散门或晶体管逻辑，离散硬件部件，或上述任何组合的设计来实现或操作所描述的 功能。通用处理器可以为微处理器，可选地，该通用处理器也可以为任何传统的处理器、 控制器、微控制器或状态机。处理器也可以通过计算装置的组合来实现，例如数字信号处 理器和微处理器，多个微处理器，一个或多个微处理器联合一个数字信号处理器核，或任 何其它类似的配置来实现。

本发明实施例中所描述的方法或算法的步骤可以直接嵌入硬件、处理器执行的软件模 块、或者这两者的结合。软件模块可以存储于RAM存储器、闪存、ROM存储器、EPROM 存储器、EEPROM存储器、寄存器、硬盘、可移动磁盘、CD-ROM或本领域中其它任意 形式的存储媒介中。示例性地，存储媒介可以与处理器连接，以使得处理器可以从存储媒 介中读取信息，并可以向存储媒介存写信息。可选地，存储媒介还可以集成到处理器中。 处理器和存储媒介可以设置于ASIC中，ASIC可以设置于用户终端中。可选地，处理器和 存储媒介也可以设置于用户终端中的不同的部件中。

在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、 固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒 介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑 存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是 任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不 限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置， 或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或 特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读 媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光 纤电缆、双绞线、数字用户线（DSL）或以例如红外、无线和微波等无线方式传输的也被 包含在所定义的电脑可读媒介中。所述的碟片（disk）和磁盘（disc）包括压缩磁盘、镭射 盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光 学复制数据。上述的组合也可以包含在电脑可读媒介中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细 说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的 保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包 含在本发明的保护范围之内。