海洋是一个时变、空变的声传输系统，海洋环境的复杂性会显著影响水声传播计算的准确性。通常海洋声学环境是未知的，并且充满不确定性，声场预报受到海洋环境不确定性效应的显著影响。声场预报中，水声传播损失的准确预报对于准确解算声呐方程，得到声呐作用距离和预测声呐性能起到关键作用。

Gerstoft等^[1]利用极大似然估计的方法尝试量化由海洋环境不确定性所导致的水声传播损失，得到水声传播损失在距离和深度两个方向上的全概率分布。Huang等^[2]进一步提出采用马尔科夫链蒙特卡洛方法利用抽样的样本将环境参数的不确定性传递给水声传播损失的计算过程，利用统计学中置信区间和相关系数对传播损失的不确定性进行量化和表征。Hoh等^[3]利用北约水下研究中心(NATO Undersea Research Center)在地中海进行的拖曳阵实验MAPEX2000数据量化由噪声和数据建模误差所引起的地声参数不确定性，并在此基础之上进行水声传播损失的统计学分析和估计。Huang等^[2]及Gerstoft等^[4]研究匹配场反演结果的不确定性对声传播损失预报的影响，通过传播损失的后验预报概率来评价传播损失预报不确定性。过武宏等^[5]通过贝叶斯理论建立水声环境不确定性推理模型，对传播损失的不确定性进行了估计，得到了传播损失80%的可信区间。也有一些学者从水声环境不确定性传播的角度去研究不确定环境参数对水声传播损失的估计的影响^[6-8]。这些方法通常在计算上比较复杂，难以在精度和效率之间取得平衡，迄今为止只有被用于简单的问题。

通常水声传播损失的后验概率密度的积分并没有解析解存在，必须采用蒙特卡洛类数值积分方法求解。在利用蒙特卡洛进行不确定性分析的过程中要反复调用声场计算模型。经典水声计算模型在水声传播计算的精度和效率上存在矛盾。直接利用传统水声模型对水声传播损失的不确定性进行分析代价过高。作为一种最优、线性、无偏估计，Kriging替代模型可以有效模拟传统水声模型输入输出响应关系，以较小的计算代价得到和模拟模型相近的输入输出关系，以平衡计算精度和效率。本文采用Kriging方法建立针对Kraken模型的声场替代模型，并集成进入蒙特卡洛不确定性分析框架，对水声传播损失的不确定性进行高效计算，在保证计算精度的同时充分提高计算效率。

1 集成Kriging模型的蒙特卡洛方法

Kriging模型最早来源于地统计学，是包含自相关(即测量点之间的统计关系)的统计模型^[9]，是一种最优线性无偏估计，可以通过平均误差的无偏差和最小化来导出。这里“最优”是指Kriging预测模型具有最小均方误差; “线性”表示Kriging预测模型可以在线性混合模型中实现，用于估计随机效应; “无偏”是指Kriging预测值的预期等于某些假设下的真实值，如高斯过程的正态性。整个模型运行需要经历一个两步过程:1)创建变异函数和协方差函数，以估算取决于自相关模型(拟合模型)的统计相关性(称为空间自相关)值；2)预测未知值(进行预测)。拟合模型或空间建模也称为结构分析或变异分析，拟合模型基于半变异函数和协方差函数。半变异函数来源于地统计分析，函数为^[10]

$ r\left( {x,h} \right) = \frac{1}{2}{\rm{Var}}\left[ {Z\left( x \right) - Z\left( {x + h} \right)} \right] $

(1)

式中:Z(x)是测量点X处的测量值；Z(x+h)是测量点x+h处的测量值。为满足平稳假设:

$ r\left( {x,h} \right) = \frac{1}{2}{\rm{E}}{\left[ {Z\left( x \right) - Z\left( {x + h} \right)} \right]^2} $

(2)

协方差函数定义为

$ {\rm{Cov}}\left( {X,Y} \right) = {\rm{E}}\left\{ {\left[ {x - E\left( Y \right)} \right]\left[ {Y - {\rm{E}}\left( Y \right)} \right]} \right\} $

(3)

在Kriging方法中，提供了多种拟合半变异模型的函数，分别是指数形式函数(EXP)、线性函数(LIN)、高斯函数(GAUSS)、三次样条函数(CUBIC SPLINE)等不同形式，针对不同应用采用不同的拟合函数，效果也不尽相同。

在对水声模型的替代过程中，Kriging模型的声场输出P由以海洋环境参数展开的多项式F和随机误差分布Z两部分组成^[10]:

$ P\left( \varepsilon \right) = F\left( {\beta ,\varepsilon } \right) + Z\left( \varepsilon \right) = f\left( \varepsilon \right)\beta + Z\left( \varepsilon \right) $

(4)

式中:F(β, ε)是回归模型，用于描述声场输出的全局近似；β为回归系数；f(ε)为含ε的多项式，可以是0阶、一阶、二阶多项式；Z(ε)为逼近误差，反映声场输出的随机性。该物理量具有如下统计特性:

$ {\rm{E}}\left( {Z\left( \varepsilon \right)} \right) = 0 $

(5)

$ {\rm{Var}}\left( {Z\left( \varepsilon \right)} \right) = \sigma _z^2 $

(6)

$ {\rm{Cov}}\left[ {Z\left( {{\varepsilon _i}} \right),Z\left( {{\varepsilon _j}} \right)} \right] = \sigma _z^2\left[ {{R_{ij}}\left( {\theta ,{\varepsilon _i}{\varepsilon _j}} \right)} \right] $

(7)

式中R_ij(θ, ε_i, ε_j)为相关函数，用以拟合建立半变异函数模型。经过对比，本文利用GAUSS相关函数和二阶多项式建立半变异函数模型。

蒙特卡罗方法又称为统计实验方法^[11]，采用蒙特卡罗方法可进行大量的统计试验，通过统计实验的结果，获得问题的近似解和概率度量，可用于不确定性分析。然而，利用经典蒙特卡洛方法进行传播损失不确定性分析时，需要反复调用声场计算模型，计算代价巨大。比如，利用Kraken模型进行声场计算，将两个地声参数作为不确定环境变量，每个参数进行1 000次蒙特卡洛模拟，则一共需要调用100万次Kraken声场计算，计算代价极高。本文利用Kriging模型建立Kraken声场计算模型的替代，并将其集成到蒙特卡洛仿真计算框架之中进行水声传播损失的不确定性量化，从而极大的减小计算代价，获得计算精度与效率之间的平衡。

2 数值计算与分析 2.1 基于Kriging的声场替代模型

Kraken是基于简正波理论求解波动方程的声传播计算程序^[12]，采用有限差分方法求解方程。利用该方法，可以快速得到问题的精确解，并可以将离散化模式问题转变为求解代数特征值问题，是水声传播计算的经典方法。本研究构建如图 1所示^[13]的环境模型并利用Kraken进行水声传播的标准计算。在实际实现中，已有学者编写了基于Matlab语言的声学工具箱。Kraken模型是该工具箱的一部分，是用于实现二维声场信道的计算模块。

实验采用两层海底模型，分别为沉积层和海底。声源位置、声速剖面等环境参数如图 1所示。本实验选取海底的密度与声速两个变量作为不确定环境变量，分别设置取值区间和分布，通过蒙特卡洛模拟进行不确定环境参数对声传播损失的影响分析。

在水声传播计算中利用Kriging模型对Kraken模型进行替代。运行Kraken模型，得到一定环境参数下固定点的声压数据，计算得到相应水声传播损失:

$ {T_\rm{L}} = - 20\lg \left| {P/P\left( {r = 1} \right)} \right| $

(8)

式中P为Kraken计算得到的声压值。在获得声传播损失之后，选取包括边界点在内的一定数量的环境参数及对应的声传播损失值作为Kriging模型的输入和输出，对Kriging模型进行训练。由于Kriging模型在训练过程中，会产生对应的估计误差。通过在估计误差(mean square error, MSE)最大值处不断添加新的训练点，逐步提高Kriging模型的精度，并最终使Kriging模型误差收敛到可接受的范围。在实际利用Kriging模型进行替代计算时，由于计算区域的空间分辨率和训练Kriging模型的复杂度存在矛盾(作为一种全局模型，Kriging模型预测的空间范围越广则其空间上的分别率越低)，分别建立了针对全部计算区域的Kriging替代模型以及以0.5 km的水平距离为分段间隔的一系列Kriging替代模型，来分别对不确定环境参数对水声传播损失的影响进行研究。全部计算区域的Kriging替代模型的环境参数设置如下:水平距离2~10 km，深度200 m，海底声速1 785~1 805 m/s，海底密度1~2 g/cm³，计算次数550次。

在对全区域建立Kriging替代模型时，总体精度较高。但由于模型本身复杂度限制，拟合点的间隔较大，空间分辨率较低。因此，可以利用该模型对某个点的传播损失后验概率密度进行求解，但不适合于对整个区域的传播损失后验概率求解。图 2为针对全部计算区域的Kriging替代模型的估计误差即最小均方误差MSE的变化趋势^[14]。图 2(a)为MSE的最大值随着训练样本增大的变化，图 2(b)为MSE的总和随着训练样本增大的变化。由图可以看出，MSE最终收敛到一个极小值，显示Kriging模型达到一个较高的精度。图 3(a)展示了全部计算区域Kriging模型最终的相对误差，由图可以看出，绝大部分误差集中在0~0.1，模型精度较高。表 1展示了Kriging替代模型估计误差的值。分段替代的Kriging模型的环境参数设置如下:水平距离1~5 km; 深度87.5、92.5、97.5 m、海底声速1 785~ 1 805 m/s、海底密度1~2 g/cm³、计算次数220次。Kriging分段替代模型基于固定深度，并以0.5 km的水平距离为间隔。由于采样点多，间隔点间距较小，因此Kriging模型的空间分辨率较高可以很好对整个区域的水声传播损失后验概率进行预测。图 3(b)展示了分段Kriging模型在4.5~5 km水平距离上相对误差值。计算结果可以看出，总体计算误差较小。具体实验数据见表 1。

2.2 水声传播损失的不确定性量化

针对单点的水声传播损失后验概率密度分布进行研究。环境参数的设置以及总计算时间(包括替代模型建立与不确定分析)如下:Kriging单点替代模型的环境参数设置如下:水平距离2~10 km、深度200 m、海底声速1 800±10 m/s、海底密度2±1 g/cm³、总计算时间324.78 s。

设定不确定海底声速与海底密度在上述区间内满足均匀分布，并将水平距离R固定为2.5、3.5和4.5 km，对应的深度D为85、90和95 m。对每一个计算点随机抽取2 000组数据进行蒙特卡洛模拟。计算点处的传播损失概率密度分布如图 4所示，频数分布如图 5所示。集成了Kriging替代模型的蒙特卡洛不确定性分析，其计算时间控制在1 s以内(计算平台为Intel i7-4790 CPU，8 GB RAM，Matlab版本为R2015b)。由计算结果可以看出，当计算深度D为85和90 m时，整体传播损失值后验概率分布较为集中。与之相对，当计算深度D为95 m时，水声传播损失后验分布较为分散；在D为85、R为2.5 km时，数据集中在一个数值附近，随着距离的增加，数据出现呈现发散的趋势。在R为95 m时，预测结果向两个值靠拢; 到R为4.5 km时，概率密度分布出现两个极值。总体来说，传播损失的不确定性在固定深度，随着传播距离的增加而增加，在同一距离，不确定性随深度变化并没有明显规律。

除了单点的传播损失后验概率密度分布计算，我们还利用分段的Kriging替代模型，对整个计算区域的传播损失进行预测，参数设置及总计算时间如下:水平距离1~5 km；深度为87.5和92.5；海底声速1 795 m/s，海底密度1±0.5 g/cm³；总计算时间173.79 s。

通过固定深度，取尽可能多的采样点来提高水声传播损失后验概率密度估计的空间分辨率。在0.5 km的水平分段里，取100个点，对每个点进行2 000次蒙特卡罗模拟。因为不确定性分析基于Kriging替代模型，其计算为线性计算，每个点的计算时间控制在1 s以内(计算平台为Intel i7-4790 CPU, 8 GB RAM，Matlab版本为R2015b)，图 6为1~5 km水平距离内水声传播损失的置信区间分布，红色线条为90%置信区间上限，蓝色线条为90%置信区间下限^[15]。

3 结论

1) 在声压计算方面，基于Kriging替代模型声场计算方法的整体计算误差控制在10%以内，计算精度较高。验证了所提出的基于Kriging替代模型的水声传播计算的可行性。

2) 通过将蒙特卡洛方法和Kriging替代模型结合，可以高效快速计算水声传播损失的概率密度分布和置信区间。总计算时间相比传统方法降低了约4个量级，获得了极大的性能提升，验证了所提出的基于Kriging替代模型的水声传播损失不确定性分析方法的有效性。

Kriging替代模型在水声传播计算当中的稳定性还有待进一步提高。下一步将针对大尺度、大范围水声传播替代计算时，Kriging模型的可靠性和稳定性进行深入分析和研究。本研究预期的研究成果可以应用到声呐性能的快速准确预报等水声应用场景。