基于海洋声速剖面数据的多插值方法拟合及误差分析

Research on Fitting Degree and Error Analysis of Multi-Interpolation Method Based on Ocean Sound Velocity Profile Data

摘要: 本文基于某次海试测点CTD获得的声速剖面数据，利用拉格朗日插值、最优点插值、样条插值和kriging插值法，得到某测点的声剖曲线，与实测数据进行拟合，通过拟合曲线误差分析，得到各插值法处理海洋声剖数据的优缺点。拉格朗日插值易在跃层区过拟合，偏离实测值。最优点插值法通过增加龙格现象抑制，能有效减少振荡。样条插值法中B样条更适配浅水区，线性与三次样条可兼顾全水深，最近邻样条平滑性最差。kriging插值拟合的声速剖面曲线整体分布上表现更均匀，出现极端误差值的概率最低，但其局部精确度上弱于其他插值方法，累计误差高且均匀分布于各个深度。

Abstract: Based on the sound speed profile data obtained by CTD at a certain sea trial station, this study employs Lagrange interpolation, optimal point interpolation, spline interpolation, and Kriging interpolation to generate sound speed profile curves for the station. These curves are fitted with the measured data, and the advantages and disadvantages of each interpolation method in processing marine sound speed profile data are derived through fitting curve error analysis. Lagrange interpolation is prone to overfitting in the thermocline region, leading to deviations from the measured values. The optimal point interpolation method can effectively reduce oscillations by adding Runge phenomenon suppression. Among spline interpolation methods, B-spline is more suitable for shallow water areas, linear and cubic splines can balance the entire water depth, and nearest neighbor spline has the poorest smoothness. The sound speed profile curves fitted by Kriging interpolation show a more uniform overall distribution and the lowest probability of extreme error values, but their local accuracy is weaker than that of other interpolation methods, with high cumulative errors evenly distributed across all depths.

1. 引言

海水声速剖面直接反映声波在海洋中的传播路径、衰减规律，对水下声学设备的性能、海洋调查任务的精度及水下作战体系保障等方面均发挥重要作用。实际海洋调查中难以获取连续声速剖面数据，因此基于有限离散的观测数据，通过插值法在已知测点数据中填补缺失的声速剖面数据，实现数据的连续性和完整性，成为海洋声学领域亟待解决的关键技术问题之一[1]-[7]。

近年来，国内外诸多学者开展了相关研究。Geyer F等[6]将拉格朗日插值用于补全离散的北极弗拉姆海峡声速观测数据。罗宇等[8]为修正超短基线定位中声速不均匀导致的声线弯曲问题，先用线性插值对声速剖面数据加密，再通过二次多项式拟合声速，构建声线跟踪模型。周芃希等[6]提出结合拖曳设备离散测点数据的声速剖面重构方法，本质是优化后的最优点插值思路，即筛选三组拖曳温深传感器的关键数据点，通过构建DHNM基函数计算系数，实现深海声速剖面重构。马海波等[9]通过距离分组试验确定Kriging插值的变异函数，将其与反距离加权法等对比，得出Kriging插值在潮间带声速关联的测深数据插值中拟合精度更高的结论。Chen P等[10]指出，传统Kriging插值在水下地形和关联声速数据处理中存在过度平滑问题，通过引入分形补偿改进后，适用于与地形强相关的近海底声速数据插值场景。

本文基于某次海上观测任务中实测的标准化1 m间隔的CTD数据，利用经典声速经验公式[11] [12]，得到声速剖面图，重点是用各插值方法进行拟合，进行曲线对比和误差检验分析，优选声剖数据的插值处理方法。

2. 原理

2.1. 拉格朗日插值法原理

过两点(x₁, y₁)、(x₂, y₂)确立一条直线，

$\begin{array}{c}\alpha \left(x\right)=ax+b\end{array}$ (1)

为确定参数a、b，把两点坐标代入公式(1)

两点式：

$\begin{array}{c}{\alpha }_{1}x=y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0}\end{array}$ (2)

一次插值基函数，分别记为$l_1=\frac{x-x_1}{x_0-x_1}$ 和$l_2=\frac{x-x_0}{x_1-x_0}$

插值函数${\alpha }_1\left(x\right)$ 是这两个基函数的线性组合形式，这种形式里的组合系数就是对应点上的函数值。此时，插值函数${\alpha }_1\left(x\right)$ 称为拉格朗日(Lagrange)插值[13]-[15]。

点斜式：

$\begin{array}{c}{\alpha }_1\left(x\right)=y_0+\frac{y_1-y_0}{x_1-x_0}\left(x-x_0\right)\\ =y_0+\frac{f\left(x_1\right)-f\left(x_0\right)}{x_1-x_0} \end{array}$ (3)

由于函数$f\left(x\right)$ 在$x_i、x_j$ 处的一阶均差为

$\begin{array}{c}f{x}_i,x_j=\frac{f\left(x_i\right)-f\left(x_j\right)}{x_i-x_j}\end{array}$ (4)

因此，(2.4)公式中的$\frac{f\left(x_i\right)-f\left(x_j\right)}{x_i-x_j}$ 是$f\left(x\right)$ 在$x_1、x_0$ 处的一阶均差$f\left[x_1,x_0\right]$ 。由于均差的对称性，${\alpha }_1\left(x\right)$ 公式的正确表述方法为

$\begin{array}{c}{\alpha }_1\left(x\right)=f\left(x_0\right)+\left(x-x_0\right)f\left[x_0,x_1\right]\end{array}$ (5)

2.2. 最优点插值法原理

计算Chebyshev节点

$\begin{array}{c}{x}_k =\cos \frac{\left(2k-1\right)}{2n}\pi ,k=1,2,\dots ,n\end{array}$ (6)

其中，$n$ 是多项式的度数加一，$x_k$ 是区间$\left[-1,1\right]$ 上的第k个Chebyshev节点。

得到Chebyshev节点后，推导最优点插值的关键在于对插值多项式的误差进行分析。对于插值多项式$P_n\left(x\right)$ ，其在某个函数$f\left(x\right)$ 上的插值误差可表示为：

$\begin{array}{c}{E}_{n\left(x\right)} = f\left(x\right)- P_{n\left(x\right)}\end{array}$ (7)

当$f\left(x\right)$ 在插值区间上具有足够的连续导数时，根据插值理论，误差可以写为：

$\begin{array}{c}{E}_{n\left(x\right)} =\frac{f^{\left(n+1\right)} \left({\xi }_x\right) }{\left(n+1\right)!}{\displaystyle \prod }_n^{i=1}\left(x-x_i\right)\end{array}$ (8)

其中，${\xi }_x$ 是一个位于$x$ 和插值节点间的某个未知点，${\displaystyle \prod }_n^{i=1}x-x_i$ 是节点$\xi$ 的有关$x$ 的多项式，代表$\frac{f^{\left(n+1\right)} \left({\xi }_x\right) }{\left(n+1\right)!}$ 的最大值。当使用Chebyshev节点时，在$\left[-1，1\right]$ 这个区间内这个多项式的值是最小的，而$f^{\left(n+1\right)}\left({\xi }_x\right)$ 表示$f\left(x\right)$ 的$\left( n+1\right)$ 阶导数在某个点$\xi$ 的值[15]-[19]。

因此，如果$f\left(x\right)$ 的$\left( n+1\right)$ 阶导数在$\left[-1,1\right]$ 区间上存在且连续，使用Chebyshev节点会使插值误差的上界最小化。

2.3. 样条插值法原理

以三次样条插值为例：假设1个数据点集$\left[\left(x_1,y_1\right),\left(x_2,y_2\right),\dots ,\left(x_n,y_n\right)\right]$ ，其中$x_1 < x_2 <\dots < x_n$ 。目标是找到一个函数S (x)，在每个$\left[x_1,x_{\left(i+1\right)}\right]$ 上是一个三次多项式，并且在整个定义域上连续的基础上，一阶和二阶导数也连续[20]-[24]。

下一步，构建构造三次多项式。在每个区间$\left[x_i,x_{\left(i+1\right)}\right]$ 上，样条函数$S_{i\left(x\right)}$ ：

$\begin{array}{c}{S}_{i\left(x\right)} = a_i + b_i \left(x - x_i\right)+ c_i {\left(x - x_i\right)}^2 + d_i {\left(x - x_i\right)}^3\end{array}$ (9)

其中$a_i,b_i,c_i,d_i$ 均为系数，需要通过边界条件和连续性条件推导，而边界条件则需计算得出样条函数在每个数据点上的函数值必须等于该点的函数值，即$S_i\left(x_i\right)= y_i$ 和$S_i\left(x_{\left(i+1\right)}\right)= y_{\left(i+1\right)}$ 。

下一步，确保一阶导数和二阶导数连续和进行边界条件的选择。确保连续需要相邻多项式的一阶导数和二阶导数在数据点处必须相等，即$S_i^{\text{'}}\left(x_{\left(i+1\right)}\right) = S_{\left(i+1\right)}^{\text{'}}\left(x_{\left(i+1\right)}\right)$ 和$S_i^{\text{'}\text{'}}\left(x_{\left(i+1\right)}\right)= S_{\left(i+1\right)}^{\text{'}\text{'}}\left(x_{\left(i+1\right)}\right)$ ；选择正确恰当的边界条件，需选择合适的条件来确定曲线两端的二阶导数。常见的选择是自然样条，其要求两端的二阶导数为零，即$S_1^{\text{'}\text{'}}\left(x_1\right)= 0$ 和$S_{\left(n-1\right)}^{\text{'}\text{'}}\left(x_n\right)= 0$ 。

最后，联立上述公式解方程组，求出所有的三次多项式系数$a_i,b_i,c_i,d_i$ 。将系数代入多项式，获得三次样条插值函数的预测模型。

2.4. Kriging插值法原理

核心在于假设数据存在空间自相关性，即空间和数值的相近和相似是依照一定的函数关系同步变化的，这种函数关系可以通过设置变异函数γ (h)来量化，它描述了样本点随着距离增加而变化的半方差[25]-[27]：

$\begin{array}{c}\gamma \left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^{N\left(h\right)}{\left(Z\left(x_i\right)-Z\left(x_i+h\right)\right)}^2\end{array}$ (10)

其中，$Z\left(x\right)$ 是位置$x$ 的随机变量，$h$ 是样本点之间的距离，$N\left(h\right)$ 是具有相同距离$h$ 的样本对数。

目标是确定估计未知点$Z\left(x_0\right)$ 处的值的最优权重${\lambda }_i$ 。这些权重应当满足无偏估计和最小方差的条件。估计值$\widehat{Z}\left(x_0\right)$ 可以表示为：

$\begin{array}{c}\widehat{Z}\left(x_0\right)= {\displaystyle \sum }_n^{i=1}{\lambda }_{i}Z\left(x_i\right)\end{array}$ (11)

为得到最优估计值进行权重${\lambda }_{\text{i}}$ 的选择，一般则是通过最小化估计值的方差${\sigma }^2$ 来实现的，即：

$\begin{array}{c}{\sigma }^2 = Var[\widehat{Z}\left(x_0\right)- Z\left(x_0\right)]\end{array}$ (12)

方差可以根据权重和变异函数展开。针对于Kriging法需要无偏、最小的方差条件，分别可得：

① 依照无偏条件，能够得到权重${\lambda }_i$ ($i$ = 1, 2, 3…)满足关系式

$\begin{array}{c}{\displaystyle \sum }_{i=1}^n{\lambda }_i=1\end{array}$ (13)

② 依照无偏条件，Kriging方差${\sigma }^2$ 为最小，可以得到求解权重${\lambda }_i$ 的方程

$\begin{array}{c}{\displaystyle \sum }_{i=1}^n{\lambda }_{i}C\left(x_i,x_j\right)\mu =C\left(x_0,x_j\right)\left(j=1,2,\dots \dots ,n\right)\end{array}$ (14)

联立公式(2.13)与(2.14)，可求得待定权系数${\lambda }_i$ 的值。

3. 插值法拟合曲线对比

基于某次海上观测任务中CTD得到的1米间隔的声速剖面调测实验数据，已经过预处理，格式统一。基于1 (11.6949 N, 95.195 E)、2 (11.4005 N, 94.9003 E)、3 (11.4007 N, 95.4959 E)、4 (11.1005 N, 95.2004 E)四个坐标测点在深度0 m~881 m区间内均匀分布的声速剖面数据片段，对处于5 (11.401 N, 95.1935 E)坐标的测点声速剖面使用拉格朗日插值法、最优点插值法(在反距离平方加权插值基础上，增加了龙格现象抑制)、样条插值法、Kriging法进行拟合。其中Kriging法常见的趋势函数分零阶、一阶、二阶多项式三种，而所使用的半变异函数则包含指数函数、广义指数函数、高斯函数、线性函数、球状函数、三次样条函数六种。对应生成各拟合数据与实测数据的曲线对比图。

Figure 1. Comparison of fitting curves by lagrange interpolation

图1. 拉格朗日插值拟合曲线对比图

Figure 2. Comparison of fitting curves by optimal point interpolation

图2. 最优点插值拟合曲线对比图

如图1所示，2站、3站拉格朗日插值法对于声速剖面的拟合曲线在0 m~100 m、150 m~300 m处显著偏离实际声速值。4个站位在0 m~400 m有显著异常偏离，极值出现在200 m处，因为200 m附近是海洋声速跃层区，4个站位拉格朗日是3次多项式拟合，跃层处声速差做高次插值，导致结果在跃层区出现过拟合。

如图2所示，在最优点插值法反距离平方加权插值基础上，增加龙格现象抑制，使得拟合曲线相比拉格朗日插值更贴近实测数据，0 m~100 m和150 m~300 m处振荡减少，更平滑。

Figure 3. Comparison of Fitting Curves by Different Spline Interpolation Methods

图3. 不同样条插值法拟合曲线对比图

如图3所示，四种样条插值方法中，线性样条与三次样条拟合趋势相近且全水深偏差规律一致，最近邻样条整体拟合最不平滑，B样条在0 m~160 m拟合度最优但中深层存在分层偏差，浅水区插值优先选B样条，全水深兼顾稳定性可选线性或三次样条。

注：分别为kriging法的零阶，一阶回归模型在不同半变异函数处理情况下拟合声剖值与原始数据的对比曲线图。其中，大写字母表示使用的半变异函数，从A到F依次为指数函数、广义指数函数、高斯函数、线性函数、球状函数和三次样条函数；小写字母代表回归模型，a为零阶，b为一阶。

Figure 4. Comparison of fitting curves for sound speed profile data by kriging method

图4. kriging法拟合声速剖面数据的曲线对比图

如图4所示，不同函数搭配下的kriging插值法在处理声速剖面数据方面均有良好的拟合效果，但随着函数的变化，拟合效果也存在差异。

4. 误差检验分析

引入平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)作为误差评价指标。其中，MAE是预测值与实际观测值之间差异的平均值，对异常值不敏感，也不能反映出误差的分布情况；MRE是预测值与实际观测值之间相对差异的平均值，与MAE不同，MRE考虑了观测值的数值，能反映估计值对于观测值的准确度：RMSE是预测值与实际观测值之间差异的平方的平均值的平方根[16] [17]，可反映估计值的灵敏度与极值。三种误差指标的表达式为：

$\begin{array}{c}MAE=\frac{{{\displaystyle \sum }}_{i=1}^{n}ABS\left(y_i-\widehat{y_i}\right)}{n}\end{array}$ (15)

$\begin{array}{c}MRE=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|\frac{ABS\left(y_i-\widehat{y_i}\right)}{y_i}\right|\end{array}$ (16)

$\begin{array}{c}RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2} \end{array}$ (17)

其中，$y_i$ 和$\widehat{y_i}$ 分别是所预测站点声速剖面在深度$i$ 处的实际观测值和预测值，$n$ 为参与误差检验的总点数。对上述各插值方法生成的曲线进行误差统计分析，得：

如表1~4所示，常见插值方法中，最近邻样条插值最优(MAE 0.474, MRE 0.031%, RMSE 0.828)，线性和三次样条插值次之，拉格朗日插值误差最大(MAE 5.198, MRE 0.345%, RMSE 7.802)。Kriging插值以零阶多项式和指数函数组合最优(MAE 2.4403, MRE 10.8515%, RMSE 0.9183)。

Table 1. Error analysis of curve fitting by common interpolation methods

表1. 常见插值方法拟合曲线的误差分析

	MAE	MRE (%)	RMSE
拉格朗日插值	5.198	0.345	7.802
最优点插值	0.770	0.051	1.181
线性样条插值	0.661	0.051	0.966
三次样条插值	0.665	0.044	1.112
最近邻样条插值	0.474	0.031	0.828
B样条插值	0.587	0.039	3.446

Table 2. MAE value of fitting data set by different functions of Kriging method

表2. Kriging法不同函数设置拟合数据的MAE值

回归模型	指数函数	广义指数函数	高斯函数	线性函数	球状函数	三次样条函数
零阶多项式	2.4403	3.5019	3.5019	3.5019	3.5019	3.5019
一阶多项式	3.5019	3.5019	3.5019	3.5019	3.5019	3.5019

Table 3. MRE value of fitting data set by different functions of Kriging method (%)

表3. Kriging法不同函数设置拟合数据的MRE值(%)

回归模型	指数函数	广义指数函数	高斯函数	线性函数	球状函数	三次样条函数
零阶多项式	10.8515	10.905	10.905	10.905	10.905	10.905
一阶多项式	10.905	10.905	10.905	10.905	10.905	10.905

Table 4. RMSE value of fitting data set by different functions of Kriging method

表4. Kriging法不同函数设置拟合数据的RMSE值

回归模型	指数函数	广义指数函数	高斯函数	线性函数	球状函数	三次样条函数
零阶多项式	0.9183	0.9183	0.9212	0.9228	0.9228	0.9228
一阶多项式	0.9206	0.9206	0.9225	0.9226	0.9226	0.9226

最近邻样条、线性、三次样条、最优点插值更适配声速剖面数据的层状连续分布特征，其MAE、MRE均处于0.031%~0.051%的低误差区间，体现了声速场空间相关性的局部化特点；而拉格朗日插值因高次多项式特性易产生龙格震荡，导致误差(MAE 5.198, RMSE 7.802)显著偏高；Kriging插值依赖数据空间相关性建模，其最优组合(零阶多项式 + 指数函数)的MAE (2.4403)和MRE (10.8515%)仍远高于最优传统插值，核心原因是声速剖面的垂直分层相关性与Kriging假设的空间连续相关性不匹配。

5. 总结

海洋声速剖面数据具有跃层突变与层状连续的特征，不同插值方法的适配性存在差异：拉格朗日插值法易在跃层区产生过拟合，导致拟合曲线显著偏离实测值；最优点插值通过反距离平方加权与龙格现象抑制，可有效降低振荡，拟合更贴近实测；样条插值中，B样条拟合浅水区更适配；线性与三次样条能兼顾全水深，更稳定；最近邻样条平滑性最差；kriging插值拟合的声速剖面曲线整体分布上表现更均匀，出现极端误差值的概率最低，但其局部精确度上弱于其他插值方法，累计误差高且均匀分布于各个深度，因此关键点是要对声速数据随深度分层，以及继续调配变异系数和回归模型。总而言之，声速剖面插值需规避高次多项式过拟合风险，浅水区优先选择B样条，全水深兼顾稳定性可选线性样条、三次样条或Kriging插值，对振荡抑制有需求时可选用最优点插值。

需要补充的是，本文结论是对于深度较浅海域的声速剖面插值处理得到的，而且不存在无条件适用的最优方法。另外，参与计算剖面数据较少，结论会受到偶然性影响。后续应在不同海洋环境下进行更多深海声速剖面数据的计算整理，反复尝试，一定能够建造更加精确的插值拟合模型，进一步提高拟合声速剖面曲线的准确性和稳定性。