1.   引　言

目前有许多关于运动介质声传播的理论与计算模型. 例如, Kirby^[1]提出了一种针对稳定的非均匀运动介质, 基于半有限元求解完全控制波动方程的方法; Hassant^[2]研究了声辐射问题, 针对活塞声源, 在低马赫数流体介质中采用了有限元方法和边界有限元法. 然而, 有限元和半有限元法难以对计算结果中的声学现象进行合理的物理解释. 另一方面, 快速场只能计算固定声源、水平分层和均匀边界的环境^[3]. Elisseeff和Schmidt^[4]利用简正波求解低马赫数、稳定流的海洋波导中的声传播问题, 并且应用匹配场进行流速反演. 然而该方法采用绝热简正波, 即假设流不会引起各阶简正波之间的耦合, 因此不适用于强穿透海底或者剖面随着深度变化剧烈的环境. 此外, 抛物方程一直是近些年的热点, Ostashev等^[5]提出了不存在等效声速近似、任意马赫数和低马赫数下的运动介质超宽角度抛物方程, 并且应用于北极声传播的仿真分析^[6]. 然而抛物方程和有限元方法一样均属于纯数值计算, 对于声学现象难以给出物理解释, 并且Ostashev等^[5]提出的抛物方程是高频近似的, 这对远距离三维计算是一个很大的挑战.

射线理论产生于光学, 且广泛地应用于声学、光学、电磁学和地震学等涉及波传播的领域研究. 射线声学把声传播看成无数条垂直于等相位面的射线的传播, 这些射线与等相位面垂直, 称为声线, 代表了声传播的方向和路径. 声线经历的时间即是声传播的时间, 携带的能量即为声能量. 因此, 尽管射线理论精度受到频率的限制, 但是其描述方法相较于简正波等模型更为简洁, 结果也更直观、清晰, 并且其计算速度快, 三维适应性良好. 此外, 射线法最主要的特点是可以计算声线到达结构(时间序列)和混响, 这些是其他全波模型无法实现的. 近年来关于射线理论的研究很少. Zabotin等^[7,8]提出了利用惠更斯波前追踪技术来计算波前, 该方法直接求解程函方程, 声线轨迹成为副产物, 该方法的鲁棒性和计算效率都高于传统的波前计算方法, 但是目前未对能量进行计算. 在声能计算方面, 传统的射线声学(无论是否考虑介质运动)都尝试求解传递方程或者Blokhintsev守恒式^[9], 需要计算由相邻声线构成的声线束横截面积, 在焦散区(声能会聚的区域)其横截面积等于0, 使得此处声能无穷大, 导致结果发散. 此外, 也无法计算完美影区(没有声线到达的区域)的声场^[10,11].

波束追踪是射线声学的延续, 同样来源于光学. 与传统射线法不同的是, 波束追踪在对从声源辐射的射线完成寻迹之后, 会在每条声线的很小邻域内求解波动方程, 计算出其波束宽度, 即传统射线理论中的波束横截面积, 同时完成幅度计算并且构建波束, 最后将各声线的贡献叠加即可. 正如文献[12-14]指出, 与传统的射线法相比, 波束追踪的结果处处为单值, 克服了焦散区的奇异点问题, 从而让计算程序鲁棒性提高, 提供精度更高的结果. 此外, 该方法还消除了进入影区声能的急剧过渡. 高斯波束追踪是目前应用最多的波束追踪方法, 其令每条声线的振幅随着远离中心射线的距离按照高斯分布衰减. 经过地震学中的改进^[15,16], 高斯波束追踪广泛地应用于大气^[13,17]和水下^[12,14]声传播. 在运动介质中, Fuerkaiti等^[18]直接将该方法应用于预测飞机噪声. 这里存在的最主要的问题是他们忽略了风速场对动态射线方程的影响, 因为目前的高斯波束追踪法仅有静止介质的表达式.

因此, 本文的目的是提出一种新的高斯波束追踪方法来求解亚音速运动介质中的声传播问题. 高马赫数环境在大气声传播中非常重要, 因为大气中风速可以达到30—40 m/s, 甚至更高, 其对声场的影响毋庸置疑. 但是在海洋中, 流速一般不会超过3 m/s, 学者们往往忽略其影响, 然而仿真结果表明海流的影响不容忽略. 本文基于Ostashev等^[5]提出的高频近似任意马赫数运动介质波动方程, 参考地震声学和水声的推导方法, 充分考虑了介质运动对声线参数的影响, 给出运动介质下的动态射线方程和声压表达式. 理论结果表明, 在非均匀的风速场或流场中, 声线束内的能量不一定守恒, 这与Blokhintsev^[9]提出的能量守恒不同. 本文对大气次声远距离传播问题进行数值仿真计算, 并与两种常用N$\times$2D近似算法计算结果进行了对比, 说明三维计算的必要性; 墨西哥湾流的仿真结果说明, 复杂海区的海流对声传播的影响不能忽略.

2.   运动介质中的高斯波束模型

2.1   理论模型

考虑三维非均匀运动介质, 声速$c({{{\boldsymbol{x}}}})$和流速${{{\boldsymbol{u}}}}({{{\boldsymbol{x}}}})$均是空间坐标${{{\boldsymbol{x}}}}$的光滑函数.

第1步需求解声线的轨迹方程, 用于获得构建波束的中心射线. 程函方程和射线轨迹方程组^[19,20]可分别写为

$${{{\boldsymbol{s}}}}\cdot{{{\boldsymbol{u}}}} + sc = 1,$$

(1)

$$\dfrac{{{\rm{d}}}{{{\boldsymbol{r}}}}}{{{\rm{d}}}t} = {{{\boldsymbol{u}}}} + c\dfrac{{{{\boldsymbol{s}}}}}{s},$$

(2)

$$\dfrac{{{\rm{d}}}{{{\boldsymbol{s}}}}}{{{\rm{d}}}t} = - \dfrac{1}{c + s^{ - 1}{{{\boldsymbol{s}}}}\cdot{{{\boldsymbol{u}}}}}\nabla{c} - {{{\boldsymbol{s}}}} \times (\nabla \times {{{\boldsymbol{u}}}}) - ({\boldsymbol s}\cdot\nabla){{{\boldsymbol{u}}}},$$

(3)

其中, ${{{\boldsymbol{s}}}} = \nabla\varPsi$表示波慢度, $\varPsi$是程函, ${{{\boldsymbol{r}}}}$是声线轨迹, $t$表示传播时间.

第2步是以声线轨迹为基础, 建立射线中心坐标系. 选择以波慢度的单位向量${{{\boldsymbol{L}}}}$为主向量建立中心坐标系, 原因有二, 一是程函方程反映的是相位变化, 相速度的方向即是波慢度方向, 二是中心坐标系在TNB标架或者Frenet标架下建立, 不难发现$({{\rm{d}}}{{{\boldsymbol{r}}}}/{{\rm{d}}}t)\cdot({{\rm{d}}}^2{{{\boldsymbol{r}}}}/{{\rm{d}}}t^2)\neq0$, ${\boldsymbol{s}}\cdot({{\rm{d}}}{{{\boldsymbol{s}}}}/{{\rm{d}}}t) = 0$. 为了与后文另一坐标系区分, 称该坐标系为相射线中心坐标系$(l, m, n)$, $l$表示沿着波慢度方向的弧长, $(m, n)$则表示场点距离中心射线的垂直距离. 在TNB标架下, 主法向量和次法向量与轨迹曲线的曲率和挠率有关, 为了避免复杂计算, 参考Červený和Hron^[21]的工作, $(l, m, n)$对应的基底分别为

$${{{\boldsymbol{L}}}} = (L_1, L_2, L_3) = {{{{\boldsymbol{s}}}}}/{s},$$

(4)

$$\begin{split} {{{\boldsymbol{e}}}}_{1L} =\, &\bigg( \dfrac{L_1L_3\cos\varTheta + L_2\sin\varTheta}{\sqrt{L_1^2 + L_2^2}}, \\& \dfrac{L_2L_3\cos\varTheta - L_1\sin\varTheta}{\sqrt{L_1^2 + L_2^2}}, \\& - \sqrt{L_1^2 + L_2^2}\cos\varTheta \bigg), \end{split}$$

(5)

$$\begin{split} {{{\boldsymbol{e}}}}_{2L} =\;&\bigg( \dfrac{L_1L_3\sin\varTheta - L_2\cos\varTheta}{\sqrt{L_1^2 + L_2^2}}, \\& \dfrac{L_2L_3\sin\varTheta + L_1\cos\varTheta}{\sqrt{L_1^2 + L_2^2}}, \\& - \sqrt{L_1^2 + L_2^2}\sin\varTheta \bigg),\\[-15pt]\end{split}$$

(6)

其中$\varTheta$表示旋转角度, 其满足微分方程

$$\begin{aligned} \dfrac{{{\rm{d}}}\varTheta}{{{\rm{d}}}l} = \dfrac{L_3(L_2 \partial c_{{\rm{p}}} / \partial x - L_1\partial c_{{\rm{p}}} / \partial y)}{L_1^2 + L_2^2}, \end{aligned}$$

(7)

并且可以得到相射线中心坐标系的拉梅系数为$(h, 1, 1)$

$$\begin{array}{l} h = 1 + \dfrac{c_{{{\rm{p}}}, m}}{c_{{\rm{p}}}}\Big|_{m = n = 0}m + \dfrac{c_{{{\rm{p}}}, n}}{c_{{\rm{p}}}}\Big|_{m = n = 0}n, \end{array}$$

(8)

其中, $c_{{\rm{p}}} = c + {{{\boldsymbol{L}}}}\cdot{{{\boldsymbol{u}}}}$是相速度, $c_{{{\rm{p}}}, m} = \partial c_{{\rm{p}}}/\partial m$和$c_{{{\rm{p}}}, n} = \partial c_{{\rm{p}}}/\partial n$表示相速度在两个方向的偏导数.

第3步是求解动态射线方程组和中心射线的声压. 动态射线方程组和中心射线的声压表达式是在中心射线附近利用抛物方程求解波动方程得到的, 详细的推导过程参考附录A, 此处仅给出最后的结果:

$$\begin{array}{l} \dfrac{{{\rm{d}}}{{{\boldsymbol{P}}}}}{{{\rm{d}}}l} = - \dfrac{{\boldsymbol{CQ}}}{c_{{\rm{p}}}^2}, \end{array}$$

(9)

$$\begin{array}{l} \dfrac{{{\rm{d}}}{{{\boldsymbol{Q}}}}}{{{\rm{d}}}l} = c{\boldsymbol{TP}}, \end{array}$$

(10)

其中矩阵${{{\boldsymbol{P}}}}$表示垂直于波慢度平面的波束慢度; 矩阵${{{\boldsymbol{Q}}}}$表示沿着声线轨迹波束的扩展, $|{{{\boldsymbol{Q}}}}|$就表示波束宽度. 三维环境下, ${{{\boldsymbol{P}}}}$和${{{\boldsymbol{Q}}}}$可以分别表示为两个线性无关解:

$$\begin{array}{l} \begin{bmatrix} {{{\boldsymbol{P}}}} \\ {{{\boldsymbol{Q}}}} \end{bmatrix} = \begin{bmatrix} {\boldsymbol{p}}_1 \\ {\boldsymbol{p}}_2 \\ {\boldsymbol{q}}_1 \\ {\boldsymbol{q}}_2 \end{bmatrix} = \begin{bmatrix} \widetilde{p}_1 & \hat{p}_1 \\ \widetilde{p}_2 & \hat{p}_2 \\ \widetilde{q}_1 & \hat{q}_1 \\ \widetilde{q}_2 & \hat{q}_2 \end{bmatrix}, \end{array}$$

(11)

其中, $\widetilde{p}_1$, $\widetilde{p}_2$, $\widetilde{q}_1$和$\widetilde{q}_2$表示波束慢度和波束宽展在$m$方向的值; $\hat{p}_1$, $\hat{p}_2$, $\hat{q}_1$和$\hat{q}_2$表示波束慢度和波束宽展在$n$方向的值. ${\boldsymbol{C}}$是相声速关于$m$和$n$的二阶偏导数与这两个方向介质马赫数的复合矩阵, 矩阵${{{\boldsymbol{T}}}}$是一个耦合矩阵,

$$\begin{array}{l} {{{\boldsymbol{C}}}} = \begin{bmatrix} (1 - M_m^2)c_{{{\rm{p}}}, mm} \\ \sqrt{(1 - M_m^2)(1 - M_n^2)}c_{{{\rm{p}}}, mn}\\\sqrt{(1 - M_m^2)(1 - M_n^2)}c_{{{\rm{p}}}, mn}\\ (1 - M_n^2)c_{{{\rm{p}}}, nn} \end{bmatrix}, \end{array}$$

(12)

$$\begin{array}{l} {{{\boldsymbol{T}}}} = \begin{bmatrix} 1 \\ -\dfrac{M_mM_n}{\sqrt{(1 - M_m^2)(1 - M_n^2)}} \\ \\-\dfrac{M_mM_n}{\sqrt{(1 - M_m^2)(1 - M_n^2)}} \\ 1 \end{bmatrix}, \end{array}$$

(13)

其中$c_{{{\rm{p}}}, ij} = \partial^2 c_{{\rm{p}}}/\partial i\partial j$, $M_i = u_i/c$, $u_i = {\boldsymbol{u}}\cdot{\boldsymbol{e}}_{i}$. 上述方程组反映了声线是如何因为声源位置变化或者声线角度变化而受到扰动的.

当流速等于0时, (9)式和(10)式与文献[14]中的方程(33)一致, 此时射线的相速度、群速度等同于声速. 当流速不为0时, 相速度、群速度不再等同于声速, 此时波束慢度变化量与相速度和流速有关, 原本静止介质表达式中的声速替换为相速度; 对于波束扩展的变化, 依然只与声速有关, 只有当$M_n$和$M_m$均不为0时, 波束宽度参量${\boldsymbol{q}}_{i}$的变化量不再与波束慢度参量${\boldsymbol{p}}_i$一一对应.

沿着声线轨迹积分, 利用速度势函数与声压的关系^[5], 中心射线的声压表达式可以写为

$$\begin{split} & p_{{\rm{beam}}}(l) = \left[ \dfrac{c}{c_{{\rm{p}}}} - \dfrac{{{\rm{i}}}u_l}{2\omega}\dfrac{{{\rm{d}}}}{{{\rm{d}}}l}(\ln|{{{\boldsymbol{Q}}}}|) - \dfrac{{{\rm{i}}}u_l\sigma}{2cc_{{\rm{p}}}^2\omega} \right]\\ & \qquad \times \dfrac{A_0\exp{ \bigg[\displaystyle\int - \sigma/(2cc_{{\rm{p}}}^2){{\rm{d}}}l \bigg]}}{\sqrt{|{{{\boldsymbol{Q}}}}|}}\exp{({{\rm{i}}}\omega\varPsi)},\\[-10pt] \end{split}$$

(14)

其中$A_0$表示声线的初始幅度, $u_l = {\boldsymbol{u}}\cdot{\boldsymbol{L}}$, $\omega$表示声源角频率, $\sigma$的表达式为

$$\begin{split} \sigma = \;&u_lc_{{\rm{p}}}\dfrac{\partial c_{{\rm{p}}}}{\partial l} + u_lu_m\dfrac{\partial c_{{\rm{p}}}}{\partial m} + u_lu_n\dfrac{\partial c_{{\rm{p}}}}{\partial n} - u_lc_{{\rm{p}}}\dfrac{\partial u_l}{\partial l}\\& - u_mc_{{\rm{p}}}\dfrac{\partial u_l}{\partial m} - u_nc_{{\rm{p}}}\dfrac{\partial u_l}{\partial n}.\\[-10pt]\end{split}$$

(15)

对于无限大均匀介质中的点源, 其声压与距离的关系为^{[11, 22]}

$$p(r) = \dfrac{\sqrt{(1 - M^2\sin^2\gamma)} - M\cos\gamma}{r(1 - M^2\sin^2\gamma)(1 - M^2)}\exp{\Big( \dfrac{{{\rm{i}}}\omega r}{c_{{\rm{g}}}} \Big)},$$

(16)

其中$M$和$\gamma$分别为介质马赫数和流速与群速度${\boldsymbol{c}}_{{\rm{g}}}$的夹角. 考虑射线法是高频近似的, 忽略方程(14)中的$\omega^{ - 1}$项, 令方程(14)与方程(16)相等, 可以得到最终的声压方程为

$$\begin{split} p_{{\rm{beam}}} =& \dfrac{A_0c^2(l)}{c_{{\rm{p}}}(l)c(0)\sqrt{|{{{\boldsymbol{Q}}}}|}}\sqrt{1 - \dfrac{M_m^2M_n^2}{(1 - M_m^2)(1 - M_n^2)}}\\&\times \exp{\left(\int\dfrac{ - \sigma}{2cc_{{\rm{p}}}^2}{{\rm{d}}}l\right)}\exp{({{\rm{i}}}\omega\varPsi)},\\[-17pt] \end{split}$$

(17)

$$A_0 = [{(1 - M^2\sin^2\gamma) - M^2\cos^2\gamma}]/{(1 - M^2)}.$$

(18)

方程(15)和(17)说明当介质运动不均匀时, $\sigma$就有可能不等于零, 可见此时声能并不全部局限于同一声线束管中, 而将受到流体运动的影响, 在传播过程中不断随流体运动而向相邻声线束管扩散. 只有在流速为0时, 声线管束内的能量才是守恒的. 而Blokhintsev^[9]则认为介质是否运动, 在声线管束内其声能均是守恒的. 不过大部分环境中, 介质运动速度梯度和声速梯度均远远小于介质速度和声速, 因此整个声能扩散项$\exp{\left[\displaystyle\int - \sigma/(2 cc_{{\rm{p}}}^2){{\rm{d}}}l\right]}$的值趋于1, 其影响非常小.

第四步是在接收位置${\boldsymbol{x}}$, 对所有波束进行加权求和, 那么该位置处的总声压场为

$$\begin{array}{l} P({\boldsymbol{x}}) = \displaystyle\iint{W(\alpha_1, \beta_1)p_{{\rm{beam}}}{{\rm{d}}}\alpha_1{{\rm{d}}}\beta_1}, \end{array}$$

(19)

这里$W$为权值函数, $\alpha_1$和$\beta_1$表示声线出射时群速度的方位角与俯仰角. 对于高斯波束, 声线束的振幅随着远离中心射线的距离按照高斯分布减小, 因此取权值函数为

$$W(\alpha_1, \beta_1) = \begin{cases} \exp{[-0.5(a^2+b^2)]}, & a + b < 4, \\ 0, & {\rm{else}}, \end{cases}$$

(20)

其中

$$\begin{split} a(m_{R}, n_{R}) =& \dfrac{{\text{δ}}\alpha_1}{c(0)}\left| \hat{q}_1m_{R}\sqrt{1 - M_{m_{R}}^2}\right.\\[-5pt]&\left. - \hat{q}_2n_{R}\sqrt{1 - M_{n_{R}}^2}\right|\Big/|{\boldsymbol{Q}}|, \\[-10pt] \end{split}$$

(21)

$$\begin{split} b(m_{R}, n_{R}) =& |\cos\alpha_1|\frac{\text{δ}\beta_1}{c(0)}\left| \widetilde{q}_1m_{R}\sqrt{1 - M_{m_{R}}^2} \right. \\& \left.- \widetilde{q}_2n_{R}\sqrt{1 - M_{n_{R}}^2}\right|\Big/|{\boldsymbol{Q}}|,\\[-10pt] \end{split}$$

(22)

这里为方便判断声线是否对场点有影响, 将场点的坐标投影到波束扩展参数$\hat{q}_1, \hat{q}_2, \widetilde{q}_1, \widetilde{q}_2$上. 这些量表示了射线位置相对于出射角度的扰动, 从而定义了声线管束的横截面积, 即波束宽度. $(m_{R}, n_{R})$是以群速度${{{\boldsymbol{c}}}}_{{\rm{g}}}$为主向量建立的新坐标系, 称为群射线中心坐标系, 令${\boldsymbol{R}} = (R_1, R_2, R_3) = {\boldsymbol{c}}_{{\rm{g}}}/|c_{{\rm{g}}}|$, 可以得到$(m_{R}, n_{R})$对应的基底为

$$\begin{split} {{{\boldsymbol{e}}}}_{1R} = \; &\bigg( \dfrac{R_1R_3\cos\varTheta + R_2\sin\varTheta}{\sqrt{R_1^2 + R_2^2}}, \\& \dfrac{R_2R_3\cos\varTheta - R_1\sin\varTheta}{\sqrt{R_1^2 + R_2^2}}, \\& - \sqrt{R_1^2 + R_2^2}\cos\varTheta \bigg), \end{split}$$

(23)

$$\begin{split} {{{\boldsymbol{e}}}}_{2R} = &\bigg( \dfrac{R_1R_3\sin\varTheta - R_2\cos\varTheta}{\sqrt{R_1^2 + R_2^2}}, \\& \dfrac{R_2R_3\sin\varTheta + R_1\cos\varTheta}{\sqrt{R_1^2 + R_2^2}}, \\& - \sqrt{R_1^2 + R_2^2}\sin\varTheta \bigg), \end{split}$$

(24)

这里并未使用相射线中心坐标系来计算场点到中心射线的距离, 因为声线轨迹才是能量实际传播的路径, 权值函数的分布平面应该位于射线轨迹的法平面内, 而非波慢度的垂直平面.

海底反射系数的表达式可以参见文献[23]中的(8.1.15 c)式, 其对声线的反射做了详细的推导. 此外, 当射线穿过或在弱连续界面发生反射时, 波束参数${\boldsymbol{P}}$会发生跳变, 静止介质跳变关系式参见文献[15, 16], 其中弱连续界面定义为声速或声速梯度不连续的界面, 在运动介质中, 只需将相关公式中的声速替换成相速度即可, 这时弱连续界面定义为相速度或相速度梯度不连续的界面.

$$\begin{array}{l} {{\boldsymbol{\bar P}}} = {{\boldsymbol{P}}} + {{\boldsymbol{Q}}}\left[ {\begin{array}{*{20}{c}} { \pm {RT_1}}&{{RT_2}}\\ { \pm {RT_2}}&0 \end{array}} \right], \end{array}$$

(25)

其中

$$\begin{split} {RT_1} = &2\left( {\dfrac{1}{{c_{{\rm{p}}}^2}}\dfrac{{\partial {c_{{\rm{p}}}}}}{{\partial m}} \mp \dfrac{1}{{\bar c_{{\rm{p}}}^2}}\dfrac{{\partial {{\bar c}_{{\rm{p}}}}}}{{\partial \bar m}}} \right)\tan {\varphi _i} \\&- \left( {\dfrac{1}{{c_{{\rm{p}}}^2}}\dfrac{{\partial {c_{{\rm{p}}}}}}{{\partial l}} \mp \dfrac{1}{{\bar c_{{\rm{p}}}^2}}\dfrac{{\partial {{\bar c}_{{\rm{p}}}}}}{{\partial \bar l}}} \right)\tan^2 {\varphi _i},\\[-10pt] \end{split}$$

(26)

$$\begin{array}{l} {RT_2} = \left( {\dfrac{1}{{c_{{\rm{p}}}^2}}\dfrac{{\partial {c_{{\rm{p}}}}}}{{\partial n}} \mp \dfrac{1}{{\bar c_{{\rm{p}}}^2}}\dfrac{{\partial {{\bar c}_{{\rm{p}}}}}}{{\partial \bar n}}} \right)\tan {\varphi _i}, \end{array}$$

(27)

式中, $\pm$和$\mp$的上下符号分别对应透射波与反射波, $\varphi _i$表示入射波的掠射角, $\bar{c}_{\rm{p}}$等含上标的参数表示发生透射或者反射之后的值.

2.2   初　值

方程(2)、方程(3)、方程(7)、方程(9)、方程(10) 构成了高斯波束追踪系统, 其初始值包括声源位置$(x_0, y_0, z_0)$、声线出射的方位角和俯仰角、以及波束慢度和波束宽度. 针对出射角度, 主要有两套设定.

第一种是设定出射角度为相速度俯仰角$\alpha$和方位角$\beta$. 这种设定可以直接开始进行追踪系统的步进, 适用于已知静止介质中指向性的声源. 不过需要额外求解每条声线初始群速度的角度和间隔, 用于计算权值函数, 通过声线轨迹方程(2), 群速度俯仰角($\alpha_1$)和方位角($\beta_1$)可以写为

$$\begin{array}{l} \alpha_1 = \arcsin\left( \dfrac{c\sin\alpha}{c_{{\rm{g}}}} \right), \end{array}$$

(28)

$$\begin{array}{l} \beta_1 = \arctan\left( \dfrac{c\cos\alpha\sin\beta + u_y}{c\cos\alpha\cos\beta + u_x} \right), \end{array}$$

(29)

由于声速$c$和流速${\boldsymbol{u}}$均与$\alpha$和$\beta$无关, 因此可以很轻松得到${{\rm{d}}}\alpha_1$和${{\rm{d}}}\beta_1$的表达式.

另一种设定则是设置群速度角度$\alpha_1$和$\beta_1$. 该设置下, 需要先求出波慢度的单位向量. 通过流及其与群速度的夹角$\gamma$可以轻松求得群速度为

$$\begin{array}{l} c_{{\rm{g}}} = u\cos\gamma + \sqrt{c^2 - u^2\sin^2\gamma}, \end{array}$$

(30)

因此波慢度单位向量可以写为

$$\begin{array}{l} {\boldsymbol{L}} = \dfrac{c_{{\rm{g}}}{\boldsymbol{R}} - {\boldsymbol{u}}}{c}. \end{array}$$

(31)

矩阵${\boldsymbol{Q}}$和${\boldsymbol{P}}$的初始值为

$$\begin{bmatrix} \;\widetilde{p}_1(0) & \hat{p}_1(0) \\ \widetilde{p}_2(0) & \hat{p}_2(0) \\ \widetilde{q}_1(0) & \hat{q}_1(0) \\ \widetilde{q}_2(0) & \hat{q}_2(0) \;\end{bmatrix}= \begin{bmatrix}\; 1 & 0 \\ 0 & 1 \\ 0 & 0 \\ 0 & 0\; \end{bmatrix} .$$

(32)

3.   数值仿真

3.1   模型正确性验证

目前没有关于运动介质声场的标准模型, 因此本文采用刚性界面半无限大均匀介质中的声传播问题作为标准进行验证, 文献[24]给出了场点的速度势函数表达式, 根据速度势函数与声压的关系, 可以发现该问题的声压场本质是直达波声压与反射波(虚源)声压的叠加:

$$\begin{array}{l} p = p(x, y, z - z_{\rm{s}}) + p(x, y, z + z_{\rm{s}}), \end{array}$$

(33)

其中$z - z_{\rm{s}}$代表直达波, $z + z_{\rm{s}}$代表反射波

令声源高度为20 m, 频率为100 Hz, 介质马赫数为0.5, 定义x轴的正负半轴分别对应顺流和逆流方向. 图1给出20 m高度处声压级理论值与高斯波束(Gaussian beam tracing, GBT)计算结果在顺逆流方向的对比. 在120 m以内, 数值计算结果与理论值相差较大, 这是因为计算时设置的开角为$\pm15^{\circ}$. 120 m以内GBT结果没有反射声线仅存在直达声, 故而单调衰减. 120 m之后, 两者的结果基本一致. 图2给出整个20 m高度的声压级伪彩图, 可以看到整个x正半轴的声压级均小于负半轴, 即顺流能量小于逆流, 这与文献[11]中远离边界的均匀介质点源声压场的结论一致. 由动态射线方程(9)和(10)以及波束方程(17)和(18)可以很好地解释这一现象, 在均匀介质中矩阵${\boldsymbol{C}} = 0$, 则波束为$p_{{\rm{beam}}} = A_0/(c_{{\rm{p}}}l)$, 随着流与传播方向夹角的增大, 相速度和波束宽度变小, 声压幅度则增大.

3.2   大气远距离次声传播

大气中风的马赫数可以达到0.5以上, 是塑造大气声场的主要因素. 对于距离较远的声源, 声速与风速随高度的强烈变化会导致声场的到达结构非常复杂, 进而影响声场结构. 次声传播距离远、能量高, 是火山爆发、台风等自然灾害的观测手段之一. 此外, 还有学者提出利用火山次声来监测热层中的风^[25,26]和利用声波来估计炸弹爆炸当量^[27]. 风场对波前和声线轨的影响在文献[7, 8] 中已经进行了分析, 特别是与等效声速法的对比, 进一步说明了等效声速法的不足. 文献[6]则分析了极地二维传播损失, 不过二维情况难以直观地反映出风场引起的三维声传播效应, 并且其忽略了垂直于传播平面的风对声传播的影响.

本文计算声场的大气数据来自怀俄明大学气象学院公开的雷达探测数据^[28]. 考虑密西西比州杰克逊(纬度$32.299^{\circ}$, 经度$-90.185^{\circ}$) 2018年12月20日12时的数据. 利用声速经验公式$c = 20.05\sqrt{T(1+0.511 q)}$^[11], 其中$T$为温度, 单位为开尔文(K), $q$为比湿度, 计算得到的声速剖面与测量的风速剖面如图3所示.

目前常用的${{N}}\times{\rm{2 D}}$等效算法有两种. 第一种方法仅考虑传播平面内风速变化, 即忽略横风的作用, 其等效风速表达式为

$$\begin{array}{l} u_{{\rm{eff}}} = u_x\cos\theta + u_y\sin\theta, \end{array}$$

(34)

其中$\theta$为计算方向的角度, 不对声速剖面做任何处理. 另一种则是等效声速法, 将等效风速直接与声速剖面相加求得新的声速剖面, 然后再利用现有的静止介质中的算法进行计算, 其声速表达式为

$$\begin{array}{l} c_{{\rm{eff}}} = c + u_x\cos\theta + u_y\sin\theta. \end{array}$$

(35)

两种方法都广泛应用于二维大气声学中^[29,30], 不过学者们已经发现, 在平流层高度及其以上, 存在强风和声线具有反转点时, 等效声速法的有效性和准确性存在较严重的问题^[31,32]. 等效风速法虽然可以解决反转点变化的问题, 不过当存在较强的横风时, 声线会发生水平偏转, 同样会导致声能出现变化. 如图4给出的在正北方向5 km高度, 3种方法计算的声压级曲线. 可以看出, GBT方法结果与等效风速的结果趋势起伏一致, 但其声压级更小; 等效声速会使得传播平面内的声线轨迹出现变化, 使得远距离传播声压级曲线的趋势与另两种方法存在较大差别, 因此对于计算运动介质中的声场, 即使是水平分层环境也应该采用三维模型进行计算以保证其准确性.

声源高度5 m, 频率5 Hz, 刚性边界, 东西与南北方向声压级伪彩图如图5所示. 声压级计算式为$20\log{\left( |p({{{\boldsymbol{x}}}})|/p_{{\rm{ref}}} \right)}$, 其中$p_{{\rm{ref}}} = 20\; \text{µ}{\rm{Pa}}$为大气参考声压级, $p({{{\boldsymbol{x}}}})$表示本文高斯波束追踪计算的结果. 该环境下声场最明显的特征就是存在多个表面声道, 分别位于东南北, 且厚度也不同. 最厚的表面声道位于正北方向, 厚度约为9 km, 而东、南两侧的表面声道厚度不足百米. 这些存在表面声道的区域, 在传播方向上风速的梯度均大于声速的梯度, 因此使得声传播模式发生了改变. 对于大角度出射的声线, 同样会受到风场的影响, 不过此时传播方向上的风速梯度小于声速梯度, 不会改变声传播模式. 此外, 不同角度声压级存在明显的差别, 说明该环境下存在三维声传播效应.

为了能更直观地看出风场带来的三维效应, 图6给出9个高度水平面的声压级结果, x轴正半轴方向为东, y轴正半轴方向为北. 可以发现, 在距离无关环境下, 声压级会随着角度发生变化, 三维效应非常明显. 此外, 随着接收高度增高, 表面声道的范围越来越小, 在100 m高度时, 表面声道范围在$-140^{\circ}$—$130^{\circ}$之间; 而5 km高度时范围仅在$45^{\circ}$—$130^{\circ}$之间; 当高度到达10 km之后, 表面声道基本消失. 利用声线轨迹进行分析, 图7给出$0^{\circ}$, $90^{\circ}$, $180^\circ$, $270^\circ$方位角出射的声线束轨迹图, 红色实线表示未经过界面反射的声线, 蓝色虚线表示经过界面反射的声线, 轨迹与图5的结果相符合. 根据方程(3), 风速梯度和流速梯度共同决定了声线的偏转方向, 将方程(3)化简可以得到

$$\begin{array}{l} \dfrac{{{\rm{d}}}{{{\boldsymbol{s}}}}}{{{\rm{d}}}t} = - \dfrac{\nabla c_{{\rm{p}}}}{c_{{\rm{p}}}}. \end{array}$$

(36)

因此不同方位角和俯仰角出射声线的相速度及相速度梯度会存在较大的差别, 这是水平分层介质中三维声传播效应的主要原因. 如图7(b)和图7(c)所示, 小角度出射的声线(蓝色虚线), 东南北3个方向上, 风速梯度大于声速梯度, 使得相速度梯度大于0, 故而声线向地面弯曲, 从而形成表面声道; 大角度声线(红色实线)的相速度受到风场的影响较小, 因此相速度梯度依然小于0, 但是在顺风时相速度梯度变大, 导致轨迹曲率变小, 而逆风时相速度梯度变小, 曲率则增大. 从图7(a)可以看出, 声线的水平偏转非常明显, 不过对于水平分层介质, 声线相速度梯度仅是坐标Z的函数, 因此波慢度的水平分量不会变化, 此时声线的偏转仅由方程(2)中引入的横风大小决定, 这种偏转不会影响波束宽度等参数, 但会使得权值函数中$a$, $b$取值变化, 进而影响声线之间的干涉, 这是水平分层介质中三维效应的次要原因, 也是等效风速法与纯三维计算方法之间存在差别的主要原因.

3.3   墨西哥湾流中的声传播

海洋中存在各种各样的流动, 如洋流、涡旋和潮汐. 马赫数一般量级为$10^{-3}$, 远远小于大气中风速的马赫数, 因此学者们通常忽略海流对声传播的影响. 然而有实验表明地转流^[33,34]和潮汐^[35]会影响声传播. 针对海洋中的尺度现象, 学者们采用重构声速剖面的方法对存在内波^[36]和中尺度涡^[37,38]等区域的声传播进行分析, 然而这些只考虑了流带来的温盐变化导致的声速变化, 并未考虑流对波动方程的改变. 此外, 实际海上试验测得的声速数据已经是被海流影响后的结果. 文献[39]利用等效声速的方法已经说明了流会定性定量地改变声传播模式, 不过他们采用的简化湾流使得流速的梯度率大于声速梯度, 而实际海洋中, 大部分时候流速的梯度小于声速梯度, 但是也会影响声传播.

本小节的数据来自哥白尼网站^[40]的观测模拟数据. 纬度范围为$-72.5^{\circ}$—$-68.5^{\circ}$, 经度$35.83^{\circ}$—$39.83^{\circ}$. 流速剖面如图8所示, x方向为经度, y方向为纬度.

数据集中仅包括温盐深数据, 需要经验公式进行拟合. Coppen^[41]给出了温度、盐度和深度与声速之间的经验公式:

$$\begin{split} c =\;& c(0) + \dfrac{16.23 + 0.0253T}{10^3}D + \dfrac{0.213 - 0.01T}{10^5}D^2 \\&+ \dfrac{(0.016 + 0.0002(S-35))(S - 35)}{10^4}TD, \end{split}$$

(37)

$$\begin{split} c(0) = \;&1449.05 + 4.57T - 0.0521T^2 + 0.00023T^3 \\&+ (1.333 - 0.0126T + 0.00009T^2)(S - 35), \end{split}$$

(38)

其中T, S, D分别为温度(℃)、盐度(千分之一)和深度(${\rm{m}}$). 该公式适用的温度范围为0—35 ℃, 盐度范围0—4.5%. 利用(37)式很容易计算得到声速, 剖面如图9所示.

海底地形如图10所示, 密度为1.5 g/cm³, 声速1700 m/s, 声吸收系数为0.5 dB/λ.

图11给出声源位于(–3.6 km, –4.6 km, 50 m), 频率50 Hz, 接收深度60 m的考虑流与不考虑流的传播损失及其差别. 尽管考虑流与不考虑流两者的声场结构相似, 难以从图11(a)和图11(b)看出区别, 但还是存在定性与定量的差别. 从图11(c) 能看出17 km以内一些角度上存在较大差别, 然而由于计算开角的原因, 这里是声影区, 声线稍有变化就会有明显的差别, 因此这不是本小节讨论的重点. 本小节主要分析距离30—40 km的区域, 该区域内两者的差别有5 dB左右.

图12给出0°和180°方位角出射声线的轨迹. 从图12(a)能发现顺流方向仅存在海底反射声线形成的会聚区, 而逆流则不存在反转形成的会聚区, 这是因为逆流方向表层声速逐渐减低, 因此反转点会升高, 当海面的声速小于声源位置声速时, 上反转点消失. 顺流时声线的跨度较小, 而逆流时较大, 原因是声源处声速梯度和流速梯度均为负, 根据方程(3), 顺流时相速度梯度的绝对值更大, 曲率变大跨度变小, 逆流时流速梯度为正, 因此相速度梯度更小, 跨度变大. 跨度的改变导致会聚区位置发生变化. 图13为60 m深度顺逆流的传播损失曲线, 可以看到顺流时会聚区向靠近声源的方向移动了大约1 km. 尽管逆流不存在反转会聚区, 但是在35 km处依然存在较强的能量会聚, 原因是海底反射声线在这里集中到达, 同样流会引起该会聚区向远离声源的方向移动0.6 km, 因此在会聚区附近会出现5 dB左右的差别. 此外, 图12(b)说明流也会引起声线的水平偏转, 海底海面的反射会进一步放大流引起的水平偏转, 不过水中的水平偏转则小得多, 其影响远不及大气中横风的影响. 此外, 图13的结果表明在会聚区之前的距离, 流对传播损失的影响较小, 从会聚区开始, 传播损失曲线开始出现较大的差别, 说明在马赫数很小的运动介质中, 需要经过足够长的声程, 才能体现流对声场的影响.

4.   结　论

本文推导了亚音速运动介质下的高斯波束追踪方程组, 建立了三维亚音速运动介质下的声传播模型. 该模型充分考虑了介质运动的三维性, 能有效地计算运动介质中的声压场, 适用于介质和地形均是三维变化的环境, 并且可以利用声线轨迹等信息解释由介质运动引起的声传播变化, 如会聚区位置发生变化和声线的偏转, 有潜力为复杂介质中的声传播问题提供理论和数值计算基础.

运动介质中的声传播与传播方向和介质运动方向密切相关, 主要影响因素是传播方向内介质运动速度的梯度, 若介质速度的梯度与声速梯度大小相当甚至大于声速梯度, 则该方向上的声传播模式会发生明显变化; 当介质速度的梯度相对于声速梯度更小时, 尽管声传播模式不会发生变化, 但是也会影响声线的跨度等参数, 使得会聚区位置有所改变, 导致区域内传播损失出现较大差别. 因而在海洋等介质运动马赫数远小于1的环境中, 复杂水文、地形环境中也应考虑介质运动对声传播的影响.

${{N}}\times{\rm{2 D}}$方法忽略了介质横向速度的影响, 在横向速度较强时, 纯三维计算结果与 ${{N}}\times{\rm{2 D}}$方法结果存在显著的差别, 此时声线在横向上的偏转会非常明显, 导致接收平面内的声能下降, 使得三维计算结果小于 ${{N}}\times{\rm{2 D}}$方法结果. 需要指出的是, 根据相关理论, 当介质马赫数趋于1时, 声能最大的传播方向趋于与介质运动垂直的方向. 因此在具有较强横风的大气环境中, 即使是水平分层介质, 也应该使用三维计算方法, 使得结果更准确. 在马赫数很小的介质中, 横向速度带来的额外声线水平偏转则小得多, 即使在地形复杂的区域, 也要经过界面反射才能看出较明显的偏转, 其对声传播的影响相较于与传播方向平行的介质运动小得多.

对于声线管束内声能存在扩散的原因, 以及与Blokhintsev能量守恒式的关系需要进一步研究.

附录A

Červený和Pšenčík^[15]给出了一种推导地震声学三维点源波束方程的方法, Porter和Bucker^[12]将其推广到水声中. 这里参考其推导方法, 给出三维运动流体介质高斯波束的推导过程$.$ 从Ostashev等^[5]提出的任意马赫数点源的速度势函数亥姆霍兹方程开始:

$$\left[ \nabla^2 + \dfrac{\omega^2}{c^2} + 2{{\rm{i}}}\dfrac{\omega}{c^2}{{{\boldsymbol{u}}}}\cdot\nabla - \dfrac{({{{\boldsymbol{u}}}}\cdot\nabla)({{{\boldsymbol{u}}}}\cdot\nabla)}{c^2} \right]\varPhi = 0,\tag{A1}$$

令其解为射线的表达式为

$$\varPhi = \phi(l, m, n)\exp{({{\rm{i}}}\omega\varPsi)}, \;\varPsi = \int{\dfrac{1}{c_{{\rm{p}}}}{{\rm{d}}}l}. \tag{A2}$$

在相射线中心坐标系下, 仅考虑中心射线附近很小邻域内的变化, 其变化主要沿着$l$方向, 在$m$和$n$方向基本无变化, 即$\partial\phi/\partial m\approx0$, $\partial\phi/\partial n\approx0$, $\partial^2\phi/\partial l\partial m\approx 0$, $\partial^2\phi/\partial l\partial n\approx 0$, 再按照$\omega^{j}$整理得到

$$F_2\omega^2 + F_1\omega + F_0\omega^{0} = 0, \tag{A3}$$

$$F_2 = h\left( \dfrac{1}{c^2} - \dfrac{1}{h^2c_{{\rm{p}}}^2} - \dfrac{2u_l}{hc^2c_{{\rm{p}}}} + \dfrac{u_l^2}{h^2c^2c_{{\rm{p}}}^2} \right)\phi, \tag{A4}$$

$$\begin{split} F_1 =& \dfrac{{2{{\rm{i}}}}}{{h{c_{{\rm{p}}}}}}\dfrac{{\partial \phi }}{{\partial l}} + {{\rm{i}}}\dfrac{\partial }{{h\partial l}}\left( {\dfrac{1}{{{c_{{\rm{p}}}}}}} \right)\phi + \dfrac{{2{{\rm{i}}}{u_l}}}{{{c^2}}}\dfrac{{\partial \phi }}{{\partial l}} - \dfrac{{2{{\rm{i}}}u_l^2}}{{h{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial \phi }}{{\partial l}}\\& - \dfrac{{{{\rm{i}}}{u_l}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial {u_l}}}{{\partial l}}\phi - \dfrac{{{{\rm{i}}}u_l^2}}{{h{c^2}}}\dfrac{\partial }{{\partial l}}\left( {\dfrac{1}{{{c_{{\rm{p}}}}}}} \right)\phi \\&- \dfrac{{{{\rm{i}}}h{u_m}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{\partial }{{\partial m}}\left( {\dfrac{{{u_l}}}{h}} \right)\phi - \dfrac{{{{\rm{i}}}h{u_n}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{\partial }{{\partial n}}\left( {\dfrac{{{u_l}}}{h}} \right)\phi,\\[-15pt] \end{split}\tag{A5}$$

$$\begin{split} F_0 =\; &h\left(1 -\dfrac{u_m^2}{c^2}\right)\frac{{{\partial ^2}\phi }}{{\partial {m^2}}} + h\left(1 -\dfrac{u_n^2}{c^2}\right)\frac{{{\partial ^2}\phi }}{{\partial {n^2}}} \\&- 2h\frac{u_mu_n}{c^2}\frac{\partial^2\phi}{\partial m \partial n} + \cdots,\\[-10pt] \end{split} \tag{A6}$$

式中$c_{{\rm{p}}}$仅与$l$有关; $c$, ${\boldsymbol{u}}$均为$(l, m, n)$的函数, 故对$F_2$中的各项进行泰勒展开, 得到

$${F_2} = - \omega h{{\boldsymbol{g}}}\dfrac{{{\boldsymbol{C}}}}{{cc_{{\rm{p}}}^2}}{{{\boldsymbol{g}}}^{{\rm{T}}}}\phi,\tag{A7}$$

其中${\boldsymbol{g}} = [\nu, \eta]$, $\nu = \omega^{1/2}m$, $\eta = \omega^{1/2}n$. 再引入射线级数

$$\phi = \displaystyle\sum\limits_{k = 0}^K {\dfrac{{{\phi _k}}}{{{{\left( {{{\rm{i}}}\omega } \right)}^k}}}} , \tag{A8}$$

对结果仅保留$\omega$的一次项, 得到与频率无关的方程:

$$\begin{split} &- {{\boldsymbol{g}}}\dfrac{{{\boldsymbol{C}}}}{{cc_{{\rm{p}}}^2}}{{{\boldsymbol{g}}}^{{\rm{T}}}}{\phi _0} + \dfrac{{2{{\rm{i}}}}}{{{c_{{\rm{p}}}}}}\dfrac{{\partial {\phi _0}}}{{\partial l}} + {{\rm{i}}}\dfrac{\partial }{{\partial l}}\left( {\dfrac{1}{{{c_{{\rm{p}}}}}}} \right){\phi _0} + \dfrac{{2{{\rm{i}}}{u_l}}}{{{c^2}}}\dfrac{{\partial {\phi _0}}}{{\partial l}} \\&- 2{{\rm{i}}}\dfrac{{u_l^2}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial {\phi _0}}}{{\partial l}} - {{\rm{i}}}\dfrac{{{u_l}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial {u_l}}}{{\partial l}}{\phi _0} - {{\rm{i}}}\dfrac{{u_l^2}}{{{c^2}}}\dfrac{\partial }{{\partial l}}\left( {\dfrac{1}{{{c_{{\rm{p}}}}}}} \right){\phi _0}\\& - {{\rm{i}}}\dfrac{{{u_m}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial {u_l}}}{{\partial m}}{\phi _0} - {{\rm{i}}}\dfrac{{{u_n}}}{{{c_{{\rm{p}}}}{c^2}}}\dfrac{{\partial {u_l}}}{{\partial n}}{\phi _0} + \dfrac{{{\partial ^2}{\phi _0}}}{{\partial \nu^2 }} + \dfrac{{{\partial ^2}{\phi _0}}}{{\partial \eta^2 }}\\& - \dfrac{{2{M_m}{M_n}}}{{\sqrt {(1 - M_m^2)(1 - M_n^2)} }}\dfrac{{{\partial ^2}{\phi _0}}}{{\partial \nu \partial \eta }} = 0. \\[-15pt] \end{split}\tag{A9}$$

进一步假设方程(A9)的解为

$${\phi _0} = A\left( l \right)\exp \left( {\dfrac{{{{\rm{i}}}{{\boldsymbol{gG}}}\left( l \right){{{\boldsymbol{g}}}^{{\rm{T}}}}}}{2}} \right)\sqrt {{c_{{\rm{p}}}}\left( l \right)}, \tag{A10}$$

可以得到两个常微分方程:

$$- c{{\boldsymbol{C}}} - cc_{{\rm{p}}}^2\dfrac{{{{\rm{d}}}{{\boldsymbol{G}}}}}{{{{\rm{d}}}l}} - {c^2}c_{{\rm{p}}}^2\left( {{\boldsymbol{GTG}}} \right) = 0, \tag{A11}$$

$$2cc_{{\rm{p}}}^2\dfrac{{{{\rm{d}}}A}}{{{{\rm{d}}}l}} + \left[ {\sigma + {c^2}c_{{\rm{p}}}^2 {{{{\rm{tr}}}}({\boldsymbol{TG}})} } \right]A = 0.\tag{A12}$$

方程(A11)是关于${\boldsymbol{G}}$的Riccati型常微分方程, 可以得到射线的动态方程组为

$${{\boldsymbol{G}}} = {{\boldsymbol{P}}}{{{\boldsymbol{Q}}}^{ - 1}}, \tag{A13}$$

$$\dfrac{{{{\rm{d}}}{{\boldsymbol{Q}}}}}{{{\rm{d}}}l} = c{\boldsymbol{TP}},\tag{A14}$$

$$\dfrac{{{{\rm{d}}}{\boldsymbol{P}}}}{{{\rm{d}}}l} = - \dfrac{{\boldsymbol{CQ}}}{{c_{{\rm{p}}}^2}}.\tag{A15}$$

最后求解方程(A12)就可以得到中心射线速度势函数的表达式为

$$\phi_{{\rm{beam}}} = \dfrac{A_0}{\sqrt{|{\boldsymbol{Q}}|}}\exp\left(\int{\dfrac{ - \sigma}{2cc_{{\rm{p}}}^2}{{\rm{d}}}l}\right)\exp{\left( {{\rm{i}}}\omega\varPsi\right)}. \tag{A16}$$