1.   引　言

基于地面观测信号对爆炸声源参数进行估计是一种被动式远距离信息感知手段, 其中声波探测是对声源能量等信息进行估计的重要手段. 地表附近爆炸产生的压力波在大气中衰减为声波传播, 传播过程中大气对其各项吸收效应与其频率呈负相关, 而次声波频率较低(0—16 Hz), 在大气中传播衰减小, 是对爆炸等事件远距离监测的有效手段. 全面禁止核试验条约组织(CTBTO)为了监测全球核试验状况, 建立了预期包含60个次声台阵的国际监测系统(international monitoring system, IMS), 并将次声作为爆炸源定位与能量估计的主要手段. 现有的远距离能量估计方法为来自于理论模型搭建与实验测量数据相结合的半经验公式. 美国空军技术应用中心(Air Force Technical Applications Center, AFTAC)根据理论模型与观测数据提出比较适合低当量实验结果的能量估计半经验公式, 使用参数为声压峰值和距离^[1]. 1995年洛斯·阿拉莫斯国家实验室(Los Alamos National Laboratory, LANL)根据声波在大气中衰减特性与多次爆炸实验的数据总结, 总结出普适性更佳, 包含声压峰值、距离与大气风速参数的能量估计半经验公式^[2]. 俄罗斯岩石动力学研究所(The Institute for the Dynamics of the Geospheres, IDG)根据传播方向与风向的垂直关系分别进行修正, 对声源能量进行估计^[3]. 2002年, Kulichkov根据脉冲信号N形波动量保持理论, 利用远距离接收的N形波与U形波的冲量对声源能量进行估计^[4].

由于大气参数呈层状分布, 温度与水平风速随海拔的变化使得次声在大气中以波导形式传播. 大气参数的垂直参数分布决定声波在大气中不同传播方向的波导模式. 在LANL提出的声源能量估计公式中, 使用平流层顶水平风速作为大气参数修正因子, 该公式由于计算速度快、误差相对较小而被广泛使用, 但其修正因子对于大气参数的利用率很低, 这导致其修正范围有限且稳健性很差. 本文从大气声传播理论出发提出一种使用大气全高度垂直参数修正的声源能量估计方法, 通过大气中计算声学方法对接收声压幅值进行修正, 并针对声源中近距离范围提出适用的超压幅值衰减模型, 大大提高了使用远距离接收次声信号幅值对声源能量估计的准确性.

2.   大气全高度垂直参数修正的能量估计方法

2.1   问题提出

LANL研究了美国内华达试验场附近观测站点全年接收声压信号随日期的变化, 观察到位于声源东侧的次声台阵接收声压值在夏季显著减小, 而位于南北方向次声台阵接收声压值全年几乎不变. 通过研究水平风速随季节的变化, 研究人员认为夏季风向与声波传播方向相反导致接收声压的异常减小, 因此使用平流层顶水平风速用于修正接收点的声压:

$${{P_{{\rm{wca}}}} = P \times {{10}^{k{v_{50}}}},}$$

(1)

其中${v}_{50}$为海拔50 km处水平风速在水平声波传播方向的分量, k = 0.019为一阶修正系数, $P$为接收声压幅度, ${P}_{\mathrm{w}\mathrm{c}\mathrm{a}}$为风速修正后的声压. 对于声波远距离传播, 研究人员采用$R/{W}^{0.5}$作为比例距离对大量化学爆炸测试数据进行回归分析, 其中$R$为传播的距离, 单位为km; $W$为声源能量, 单位为kt. 通过最小二乘法得到声源能量估计公式:

$${{\rm{log}}{P_{{\rm{wca}}}} = 3.37 + 0.68\;{\rm{log}}W - 1.36\;{\rm{log}}R.}$$

(2)

(2)式使用50 km海拔处风速以简化运算. 当风速为正值时, 平流层顶部的有效声速极大值大于地面声速, 产生平流层波导^[5], 这使得接收点声压幅度偏大; 当风速为负值时, 则次声在大气中的传播只存在热层波导, 由于热层中声波传播的经典衰减与转动衰减^[6]带来的衰减量非常明显, 这使得接收点的声压幅度偏小. 在平流层顶未出现异常扰动时, (2)式可以对声源能量进行较为准确的估计. 但在实际使用中, 使用单一参数${v}_{50}$作为大气修正因子导致声源能量估计值误差较大, 而计算大气声学综合大气温度、风速和黏滞系数等因素的影响, 可降低对声源能量估计的误差.

常用的计算大气声学方法有射线法、抛物方程法、非线性渐进波动方程法(nonlinear progressive wave equation, NPE)和有限元方法^[7,8]. 射线法使用的声线模型无法读取出任意位置点的能量值, 首先予以排除. 由于本方法需要获取观测点位置声波能量损失值, 应尽量多考虑大气中的影响因素, 而有限元方法在远距离仿真状况下计算量过大, 缺乏实际可用性, 因此选用相对于抛物方程法额外考虑到折射效应与非线性陡峭等非线性效应的NPE数值模拟方法.

2.2   非线性渐进波动方程

声波在大气中的传播受大气流体特性的影响, 将声压方程组中的运动方程添加介质黏度项, 并将物态方程中的声压展开至二阶^[9]:

$$\left\{\begin{aligned} &\rho \left[\frac{{\boldsymbol{v}}}{\partial t}+\left({\boldsymbol{v}}\cdot \nabla \right){\boldsymbol{v}}\right]= \\ & -\nabla P+\left(\eta +\frac{4}{3}\mu \right)\nabla \left(\nabla \cdot {\boldsymbol{v}}\right) -\eta \nabla \times \nabla \times {\boldsymbol{v}},\\ &\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho {\boldsymbol{v}}\right)=0,\\ &P={c}^{2}{\rho }'+\frac{1}{2}{\left(\frac{{\partial }^{2}P}{\partial {\rho }^{2}}\right)}_{s}{\rho }^{{\prime}2}+{\left(\frac{\partial P}{\partial \rho }\right)}_{\rho }{s}',\end{aligned}\right.$$

(3)

其中$P$为声压, $\rho$为大气密度, s为热力学中的熵, $\eta$为体积黏度, $\mu$为剪切黏度. 相对于抛物方程等方法, NPE通过对绝热方程的展开保留二阶项以考虑声波在大气中传播的非线性效应. 将

$$\begin{split} &{\left(\frac{{\partial }^{2}P}{\partial {\rho }^{2}}\right)}_{s}=2\left(\beta -1\right)\frac{{c}^{2}}{{\rho }_{0}},\\ &{\left(\frac{\partial P}{\partial \rho }\right)}_{\rho }{s}'=\kappa \left(\frac{1}{{C}_{v}}-\frac{1}{{C}_{p}}\right)\nabla \cdot {\boldsymbol{V}},\end{split}$$

(4)

代入方程组(3)中的物态方程, 得到方程

$$\begin{split} & {\left({D}_{t}-{c}_{0}{\partial }_{x}\right)}^{2}\rho ' \\ =\;& {\partial }_{x}^{2}{\left({c}_{0}+{c}_{1}\right)}^{2} +\left({\rho }'+\beta \frac{{{\rho }'}^{2}}{{\rho }_{0}}\right)\left({\partial }_{y}^{2}+{\partial }_{z}^{2}\right) \\ & \times \Bigg\{\left({c}_{0}+{c}_{1}\right)^{2}\Bigg[\rho ' +\left(\beta -1\right)\frac{{\rho }^{{\prime}2}}{\rho_0}\Bigg]\Bigg\} \\ & -\xi {\nabla }^{2}\left[{\partial }_{x}\left({c}_{0}+{c}_{1}\right){\rho }'\right],\end{split}$$

(5)

其中热黏滞系数$\xi$为

$$\xi ={\rho }_{0}^{-1}\left[\kappa \left(\frac{1}{{C}_{v}}-\frac{1}{{C}_{p}}\right)+\left(\eta +\frac{4}{3}\mu \right)\right],$$

(6)

其中$\kappa$为热传导系数, ${C}_{v}$为定容比热, ${C}_{p}$为定压比热. 热黏滞项的计算主要估计了声波在大气中传播的经典衰减.

将方程(5)展开后, 忽略高阶小量后可化简得到NPE在笛卡尔坐标下的标准形式:

$$\begin{split} &{D}_{t}R+\frac{\partial }{\partial x}\left({c}_{1}R+\frac{{c}_{0}}{2}\beta {R}^{2}\right)-\frac{1}{2}\xi {\partial }_{x}^{2}R\\ &+\frac{{c}_{0}}{2}\int _{{x}_{\mathrm{f}}}^{x}\left({\partial }_{y}^{2}+{\partial }_{z}^{2}\right)R\mathrm{d}x=0,\end{split}$$

(7)

其中${D}_{t}={\partial }_{t}+{c}_{0}{\partial }_{x}$为运动框算子, 使得计算域随着声波传播而运动, $R={\rho }'/{\rho }_{0}$为归一化密度扰动, ${c}_{0}$为静态大气中声速, ${c}_{1}$为随运动计算框变换的声速扰动, $\beta$为非线性系数, $\xi$为热黏滞系数. 实际处理中, 考虑到大气各项参数呈层状分布, 通常使用圆柱坐标系进行数值模拟^[10]:

$$\begin{split} &{D}_{t}R+\frac{{c}_{0}}{2}\frac{R}{r}+\frac{\partial }{\partial r}\left({c}_{1}R+\frac{{c}_{0}}{2}\beta {R}^{2}\right)\\ &\frac{1}{2}\xi \left(\frac{{\partial }^{2}R}{\partial {r}^{2}}+\frac{1}{r}\frac{\partial R}{\partial r}\right)+\frac{{c}_{0}}{2}\int _{{r}_{\mathrm{f}}}^{r}\left(\frac{{\partial }^{2}R}{\partial {z}^{2}}\right)\mathrm{d}r=0.\end{split}$$

(8)

使用算子分裂法^[11]将(4)式分解为4个步骤, 并分别使用不同格式进行数值求解, 其中

$${D}_{t}R=-\frac{{c}_{0}}{2}\frac{R}{r}$$

(9)

为圆柱扩散效应项, 使用Crank-Nicolson差分格式进行数值计算;

$${D}_{t}R=-\frac{\partial }{\partial r}\left({c}_{1}R+\frac{{c}_{0}}{2}\beta {R}^{2}\right)$$

(10)

为折射与非线性效应项, 将此两项结合使用通量校正法^[12]进行数值计算可以提高计算效率;

$${D}_{t}R=\frac{1}{2}\xi \left(\frac{{\partial }^{2}R}{\partial {r}^{2}}+\frac{1}{r}\frac{\partial R}{\partial r}\right)$$

(11)

为热黏滞吸收效应项,

$${D}_{t}R=-\frac{{c}_{0}}{2}\int _{{r}_{\mathrm{f}}}^{r}\left(\frac{{\partial }^{2}R}{\partial {z}^{2}}\right)\mathrm{d}r$$

(12)

为衍射效应项, 均使用Crank-Nicolson差分格式与Thomas算法结合求解. 设定地面为刚性边界条件, 使用完美匹配层^[13]设定大气上边界为吸收边界条件, 通过运动框滑动叠加计算出声源在设定传播方向的传播损失分布.

由于NPE算法中使用的传播损失量由运动窗叠加求得, 则传播损失的绝对数值与运动框步长负相关, 无法直接用于声源能量估计, 因此需要使用传播损失的相对数值通过间接法进行计算. 根据本节前文的描述, NPE方程对非强非线性线性声波传播范围内可进行有效的模拟, 对于更近范围内的声波能量传播损失值的模拟, 由于未能保留更加高阶的小量^[14], 无法对此范围内冲击波的强非线性效应进行准确模拟. 因此, 如果直接使用接收台阵位置的传播损失量作为对声源能量估计, 声源附近的强非线性效应模拟误差与输入声源函数特性的影响为能量估计带来额外的误差, 因此本文采用弱非线性效应区域传播损失相对值法进行声源能量估计.

假设次声台阵距离声源距离为${r}_{\mathrm{t}}$, 通过NPE方程求得台阵位置传播损失量为$P{L}_{{r}_{\mathrm{t}}}$, 接收到声压单峰幅值为${p}_{{r}_{\mathrm{t}}}$. 设定距离声源${r}_{0}$为参考距离, 此位置传播损失量为${PL}_{{r}_{0}}$. 在爆炸产生的冲击波研究中通常使用的参数为比例距离$Z=r/{W}^{\tfrac{1}{3}}$, 其中r为水平距离, 单位为m; $W$为声源能量, 单位为kg^[15]. 当与声源距离满足$Z < 20 \; \mathrm{m}/{(\mathrm{k}\mathrm{g}}^{\tfrac{1}{3}})$时, 该范围为近场区域. 在此区域内爆炸声源产生的高温、高压燃爆物急剧膨胀, 将附近空气迅速排挤形成压缩空气层. 此压缩空气层以超音速运动且其状态参数有突跃, 称为冲击波, 其前沿称为波阵面. 在近场区域时, 冲击波超压幅度满足的方程满足强非线性条件. 随着传播距离的增加, 波阵面上的压强等参数快速下降, 接近衰减为声波时满足弱非线性条件. 由于本方法使用NPE作为数值模拟方法, 则参考距离应超出强非线性影响范围, 即${Z}_{0}={r}_{0}/{W}^{\tfrac{1}{3}} >$$20\; \mathrm{m}/({\mathrm{k}\mathrm{g}}^{\tfrac{1}{3}})$^[16]. 将$r$处声压转化为${r}_{0}$处超压幅值${P}_{{r}_{0}}$的转移方程为

$${P}_{{r}_{0}}={p}_{{r}_{\mathrm{t}}}\cdot {10}^{\tfrac{P{L}_{r}-P{L}_{{r}_{0}}}{20}}.$$

(13)

使用通过转移方程(13)计算得到的参考超压幅值对声源能量进行进一步估计.

2.3   中近程超压幅度衰减模型

使用2.2节计算得到的参考声压, 可根据以下超压幅度衰减模型对声源能量进行估计^[15,17,18]:

$${P}_{{r}_{0}}=\frac{8.19\times {10}^{7}\cdot \bigg[1+{\left(\dfrac{Z}{4.5}\right)}^{2}\bigg]}{\sqrt{1+{\left(\dfrac{Z}{0.048}\right)}^{2}}\sqrt{1+{\left(\dfrac{Z}{0.32}\right)}^{2}}\sqrt{1+{\left(\dfrac{Z}{1.35}\right)}^{2}}},$$

(14)

$${P}_{{r}_{0}}=1\times {10}^{6}\cdot \left(\frac{0.076}{Z}+\frac{0.255}{{Z}^{2}}+\frac{0.65}{{Z}^{3}}\right),$$

(15)

$${P}_{{r}_{0}}=1\times {10}^{6}\cdot \left(\frac{0.082}{Z}+\frac{0.265}{{Z}^{2}}+\frac{0.686}{{Z}^{3}}\right),$$

(16)

其中${P}_{{r}_{0}}$为${r}_{0}$处超压幅值. 上述超压值随比例距离变化的公式主要适用于比例距离较小的状况, 根据2016年Kim和Rodgers^[14]适用的实验测量, (14)—(16)式在$20\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}}) < Z < 50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$范围内与实测数据符合状况很好, 但在$Z > 50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$范围内, 预测值与实际测量结果偏差很大. 为了验证此类声源能量估计公式对于不同比例距离的适用状况, 本文使用Humming Roadrunner (HRR)实验数据与美国1951—1958年在内华达试验场(Nevada test site, NTS)附近各台阵测得的直达波信号幅度数据对其进行评估测试^[19], 对比结果见图1. 根据对比结果, 在$Z < 50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$范围内, 传统超压幅值衰减模型与实际测量信号符合状况较好, 反映出冲击波超压幅值随比例距离增加的快速衰减过程; 但是在$Z > 50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$范围, 实际观测到的信号随距离衰减速度明显快于传统衰减模型. 在此范围内, 传统超压幅度衰减与$1/Z$近似呈线性相关, 观测信号超压幅度亦呈线性相关, 但其系数取值差异很大.

本文提出的方法中, 参考距离${r}_{0}$应位于$Z >$$50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$范围以减少数值模拟带来的误差, 因此需要针对此比例距离范围, 根据实测数据提出合理的超压幅度随比例距离衰减的经验公式. 通过最小二乘法在对数坐标系下对实测数据进行线性拟合, 得到适用于本方法的中近程超压幅度衰减模型(middle range model, MR):

$${P}_{{r}_{0}}=4.5\times {10}^{5}\cdot {Z}^{-1.54},$$

(17)

化简为声源能量估计方程为

$$W=9.9185\times {10}^{-12}\cdot {r}_{0}^{3}\cdot {P}_{{r}_{0}}^{1.94}.$$

(18)

从图1中的NTS观测数据可以看出, 随着比例距离Z的增加, 对于超压幅值预测值的方差逐渐增加, 为了获得更高的声源能量估计精度, ${r}_{0}$取值应尽量使得对应${Z}_{0}$取值接近$50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$. 但在实际使用过程中, 本文使用NPE方程进行数值模拟时在水平方向的步进为100 m, 而NPE方法对于初始声源函数附近模拟值收敛速度较慢产生计算误差, 当${r}_{0}$取值过小时会导致声源能量估计错误. 根据多次实验结果, ${r}_{0}$应大于50倍NPE中水平方向步进, 即${r}_{0} > 5 \; \mathrm{k}\mathrm{m}$, 综合${Z}_{0} > 50\; \mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$的限制条件, 本方法仅适用于声源能量$W < 1~\mathrm{k}\mathrm{t}$的状况. 图1给出了(14)—(16)式分别对应的超压幅值衰减模型与本文提出的衰减模型效果对比, 其中KG85模型为Kinney与Graham于1985提出的化学爆炸超压衰减模型^[15]; Sadovskyi模型为根据模型相似律建立的公式, 由实验测量确定系数, 得到高爆炸药冲击波峰值超压的变化模型, 使用范围为比例距离Z满足1 < Z < 15^[17]; 王儒策模型为根据原子爆炸的经验, 球形装药在无限空气介质中爆炸时的超压变化规律^[18].

3.   实验与讨论

为了验证本计算方法的普适性, 使用2011年1月24日在以色列Sayarim进行的化学爆炸实验数据^[20]进行验证, 声源位置坐标为34.81668°E, 29.99555°N, 实验测量得到的数据与相关参数见表1, 其中的风速数据来自HWM14大气风场模型. 实验声源为在地表进行的化学爆炸, 其能量为7.4吨TNT. 由于并未获得该数据的波形数据, Kulichkov^[4]提出的基于脉冲声冲量保持理论的计算方法无法使用, 故本文将使用LANL提出的普适性能量估计公式进行对比计算, 其显式声源能量估计形式为

$$W=2.16\times {10}^{-6}\cdot {r}^{2}\cdot {p}_{r}^{1.42}{\cdot 10}^{-0.0256{v}_{50}}.$$

(19)

而在使用本文提出的基于大气传播模拟的声源能量估计方法时, 采用${r}_{0}=7\; \mathrm{k}\mathrm{m}$作为参考距离, 以HWM14与MSISE00模型作为大气参数剖面模型, 得出如图2的声源向各台阵方向的传播损失分布图.

为了验证参考距离的选取对声源估计精度的影响, 选取${r}_{0}=\mathrm{1, 2}, 3, \cdots, 30\; \mathrm{k}\mathrm{m}$进行对比计算, 误差对比见图3. 图中的变化曲线表明, 随着参考距离的减小, 除JO_NE和KU点位数据的计算结果外, 对声源能量的估计误差逐渐减小, 并在参考距离5—10 km之间达到最小误差. 对于JO_NE和KU点位数据, 由于参考声压系统性偏小, 因此计算得到的声源能量估计值的误差随着参考距离的减小逐渐增加.

选取参考距离分别为7与20 km进行声源能量估计数值的具体对比, 结果见表2, 表中${r}_{0}= 7\; \mathrm{k}\mathrm{m}$下声源能量估计结果的误差远小于${r}_{0}= 20\; \mathrm{k}\mathrm{m}$下的误差, 对应比例距离分别为${Z}_{0}=359\;\mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$和${Z}_{0}=1026\;\mathrm{m}/(\mathrm{k}{\mathrm{g}}^{\tfrac{1}{3}})$, 此误差变化结果与图1中方差变化趋势一致. 对于JO_NE台阵与KU台阵的数值模拟结果, 其中JO_NE在433 km附近传播损失很小, 此项模拟偏小导致参考声压的计算值偏小, 影响计算结果; KU台阵在此距离传播损失很小, 此项偏小导致参考声压的计算值偏小, 导致计算结果整体偏小. 使用不同的格点长度计算后, 模拟的结果近似, 表明该误差并非由格点精度导致.

根据图2中的数值计算结果, 在方位角为14°, 50°, 95.5°和113°这4个方向的次声传播过程中均存在平流层波导. 与此一致的是, 在表1中这6个点位对应的${v}_{50}$均为正值, 即认为存在平流层波导, LANL公式对声源能量估计值修正减小. 对该数据分别使用LANL提出的能量估计公式和本文提出的NPE_MR方法对声源能量进行估计, 估计结果与误差的对比结果见表3.

表3中的6个观测点位分布在距离声源300—2500 km之间的范围内, 通过表内误差的对比可以观察到, 对于此次化学爆炸试验, LANL声源声量估计公式对声源能量的估计明显偏大. 表1中的${v}_{50}$数值表明该公式对声源能量的初始估计值的修正方向正确, 但是由于大气中该高度处的风速幅值有限, 对于此数据对应的状况无法进行正确的修正. 而本文提出的NPE_MR声源能量估计方法充分利用了声波传播过程中的大气剖面参数, 对于各接收点的信号幅度进行了合适的修正, 大大降低了对爆炸声源的能量估计误差.

在实际使用中, 次声信号可以被相对于声源位置不同方位角的多个次声台阵接收到, 而LANL提出的声源能量估计公式在计算比较快速的情况下, 对于声源能量的估计误差一般在300%以内, 可直接使用多阵平均值对声源能量进行估计. 但是对于大气参数对接收信号影响较大的情况下, 传统声源能量估计方法的估计误差可能达到将近4000%, 在多阵计算中只能作为错误项去除. 而通过本文提出的NPE_MR声源能量估计方法, 可以将大气参数对接收声波幅度的影响进行大幅度修正, 得到可用于多阵平均计算的结果, 提高接收阵的信息利用效率, 减小对声源能量估计的误差.

4.   结　论

对于通过大气远距离传播的次声信号, 美国与俄罗斯使用大气中的风场信息对其进行修正以减小声源能量估计的误差. 其中美国LANL提出的包含风速修正项的声源能量估计公式对声源能量的计算误差最小, 但是在传播损失特别大时, 由于修正项可提供的修正幅度有限, 无法将估计值误差修至300%以下. 本文提出了利用大气中传播模拟仿真与MR衰减模型的声源能量估计方法, 大大拓宽了大气参数对于声源能量估计值的修正范围. 在实验验证中, 本文使用2011年以色列Sayarim化学爆炸实验的实测数据, 对比了本方法与传统半经验计算公式计算结果之间的误差, 证明本方法在需要使用大气数据对初始声源能量估计结果进行大幅度修正的情况下比传统半经验公式有更好的性能. 在后续的工作中, 需要对NPE方法的稳定性做进一步的提高, 减少声源能量估计值异常偏小的情况.