为了快速获取大面积海域的声速剖面信息，声学反演方法得到了国内外学者的广泛关注。其中匹配场反演^[1]是反演浅海声速剖面的常用方法，然而匹配场技术对环境参数的失配敏感，使得一些学者致力于匹配场反演的环境参数失配问题^[2-4]。而后张忠兵等^[5]利用波束匹配的方法与常规匹配场反演进行比较，发现波束匹配反演对底质参数失配具有较高的鲁棒性。另一类较常用的方法则是利用声线传播时间反演海水的声速剖面^[6]。张维等^[7-10]提出了基于微扰法的快速声速剖面反演，在此基础上分析了基阵倾斜失配的影响，并利用垂直阵不同阵元接收信号到达时间实现了三维声速剖面反演。由于微扰法可以实现实时的声速剖面反演，因此也被用于海流的反演中^[11-20]。

当利用声波传播时间反演声速剖面时，由于海上实验过程中信号传播时间的测量误差是难免的，并且小的测量误差有时会给声速剖面的反演精度带来极大影响。为了降低反演算法对传播时间测量误差的敏感性，本文在微扰法反演声速剖面的基础上，通过调整平均声速剖面以增大声速扰动量，以此增大平均声速剖面和实际声速剖面下的声线传播时间差，降低传播时间测量误差的影响，提高反演算法对声速扰动的敏感性，进而提高声速剖面反演算法的精度。文中对改进微扰近似所带来的反演声速剖面偏移进行修正，有效降低了平均声速剖面调整引入的误差，进一步提高了反演结果的精度。

1 微扰法改进方法

浅海环境下，经验正交函数 (EOF) 是描述浅海声速剖面最有效的基函数。当已知平均声速c(z) 和各阶经验正交函数f_k(z) 后，实际声速剖面可以用下式表示：

$ c\left( z \right) = \bar c\left( z \right) + \Delta c\left( z \right) = \bar c\left( z \right) + \sum\limits_{k = 1}^K {{\alpha _k}{\mathit{\boldsymbol{f}}_k}\left( z \right)} $

(1)

此时声速剖面反演转变为经验正交函数系数α_k的求解。微扰法计算声线的传播时间如下

$ \begin{array}{*{20}{c}} {t = \int\limits_l {\frac{1}{{\left( {\bar c\left( z \right) + \Delta c\left( z \right)} \right)}}{\rm{d}}s} \approx \int\limits_l {\frac{{\bar c\left( z \right) + \Delta c\left( z \right)}}{{{{\bar c}^2}\left( z \right)}}{\rm{d}}s} = }\\ {\int\limits_l {\frac{1}{{\bar c\left( z \right)}}{\rm{d}}s} - \int\limits_l {\frac{{\Delta c\left( z \right)}}{{\bar c\left( z \right)}}{\rm{d}}s} = \int\limits_l {\frac{1}{{\bar c\left( z \right)}}{\rm{d}}s} - }\\ {\int\limits_l {\frac{{\sum\limits_{k = 1}^{n{\rm{Vectors}}} {{\alpha _k}{\mathit{\boldsymbol{f}}_k}\left( z \right)} }}{{{{\bar c}^2}\left( z \right)}}{\rm{d}}s} = \bar t - \Delta t} \end{array} $

(2)

式中：l表示声线轨迹。当声速扰动量足够小时，l可用平均声速下的声线轨迹l₀代替。此即为微扰近似的条件：

$ \Delta t = \sum\limits_{k = 1}^{n{\rm{Vectors}}} {{\alpha _k}\int\limits_l {\frac{{{\mathit{f}_k}\left( z \right)}}{{{{\bar c}^2}\left( z \right)}}{\rm{d}}s} } = \sum\limits_{k = 1}^{n{\rm{Vectors}}} {{\alpha _k}{d_k}} $

(3)

考虑n个阵元的水平阵，声波到达各阵元的传播时间用t_i(i=1, 2, …, n) 表示，则平均声速剖面下的声线传播时间表示为

$ {t_i} = {{\bar t}_i} - \sum\limits_{k = 1}^{n{\rm{Vectors}}} {{\alpha _k}{d_{ki}}} \;\;\;\;\;\left( {i = 1,2, \cdots ,n} \right) $

(4)

写成矩阵形式：

$ \mathit{\boldsymbol{DK}} \cdot \mathit{\boldsymbol{A = \bar T}} - \mathit{\boldsymbol{T = }}\Delta \mathit{\boldsymbol{T}} $

(5)

上述由n个方程组成的线性方程组含有k个未知数。当n>k时，方程组则为超定方程组，用伪逆法求解可得

$ \mathit{\boldsymbol{A = }}{\left( {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{DK}}} \right)^{ - 1}}\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\Delta \mathit{\boldsymbol{T}} $

(6)

式 (6) 即为微扰法反演声速剖面的求解公式。由于实际声速剖面c(z) 和平均声速剖面c(z) 相近，因此平均声速下声信号传播时间t与实际声速下的传播时间t也非常接近，声速扰动所引起的传播时间扰动Δt甚小。理论计算发现，浅海中当声波水平传播距离在几十公里时，Δt一般仅为10^-3 s的量级。即便由于爆炸声源爆炸位置不准确、基阵倾斜等原因导致实测声信号传播时间仅有10^-4 s量级的微小误差，则相对Δt的误差也将达到十分之一，这对经验正交函数系数α_k的反演的影响是显著的。因此当Δt较小时，其受声信号传播时间测量误差的影响是巨大的，对声速剖面的反演是不利的。

为此，本文考虑通过改变实际声速剖面的表示形式，即调整平均声速剖面来改变平均声速与实际声速之差，以增大Δt，从而提高Δt的抗干扰能力，最终提高反演算法的鲁棒性。

引入非负常数a，调整后的平均声速表示为c′(z)=c(z)-a，这样便增大了反演中的平均声速与实际声速之差，也即增大了Δt。以c′(z) 作为新的平均声速，重新定义声速剖面数据集中的每条声速剖面与平均声速剖面之差：

$ \Delta {{c'}_i}\left( z \right) = {c_i}\left( z \right) - \bar c'\left( z \right) $

(7)

式中：c_i(z) 表示声速剖面数据集中第i条声速剖面。

将声速剖面差Δc′_i(z) 等间隔内插为N个深度，则Δc′_i(z) 离散化后写为

$ \Delta {{c'}_i}\left( {{z_j}} \right) = {c_j}\left( {{z_j}} \right) - \bar c'\left( {{z_j}} \right)\;\;\left( {j = 1,2, \cdots ,N} \right) $

(8)

类似于协方差矩阵的定义，根据离散化的M条声速剖面和平均声速剖面之差Δc′_i(z) 得到新矩阵：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{R' = }}\frac{1}{M}\left[ {\begin{array}{*{20}{c}} {\Delta {{c'}_1}\left( {{z_1}} \right)} & \cdots & {\Delta {{c'}_M}\left( {{z_1}} \right)}\\ \vdots & {} & \vdots \\ {\Delta {{c'}_1}\left( {{z_N}} \right)} & \cdots & {\Delta {{c'}_M}\left( {{z_N}} \right)} \end{array}} \right].}\\ {\left[ {\begin{array}{*{20}{c}} {\Delta {{c'}_1}\left( {{z_1}} \right)} & \cdots & {\Delta {{c'}_M}\left( {{z_N}} \right)}\\ \vdots & {} & \vdots \\ {\Delta {{c'}_i}\left( {{z_N}} \right)} & \cdots & {\Delta {{c'}_M}\left( {{z_N}} \right)} \end{array}} \right]} \end{array} $

(9)

将此矩阵R′进行特征值分解，则可以用N个正交向量表示R′，有

$ \mathit{\boldsymbol{R'F' = A'F'}} $

(10)

$ \mathit{\boldsymbol{R'}} = \sum\limits_{n = 1}^N {{{\lambda '}_n}{{f'}_n}f_n^{'H}} $

(11)

$ \mathit{\boldsymbol{A'}} = \left[ {\begin{array}{*{20}{c}} {{{\lambda '}_1}} & 0 & \cdots & 0\\ 0 & {{{\lambda '}_2}} & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots & 0 & {{{\lambda '}_N}} \end{array}} \right] $

(12)

式中：λ′_n为协方差矩阵R′的特征值，f′_n为特征值λ′_n对应的特征向量。

将协方差矩阵R′的N个特征值λ′_n按从大到小的顺序排列，即λ′₁>λ′₂>…>λ′_N，选取前K个特征值对应的特征向量用以表示声速剖面，则声速剖面可以用K阶特征向量 (广义EOF) 表示为

$ c\left( z \right) = {{c'}_i}\left( z \right) + \sum\limits_{n = 1}^N {{{\alpha '}_k}{{\mathit{\boldsymbol{f'}}}_k}\left( z \right)} $

(13)

式中：c′(z) 为调整后的平均声速剖面，α′_k为经验正交函数系数，f′_k(z) 为广义EOF。以式 (13) 作为新的声速剖面的经验正交函数表示，并利用微扰法反演声速剖面。张维^[12]针对反演精度的需求对选取不同阶的经验正交函数进行了分析，一般情况下3~4阶的经验正交函数即可满足精度要求。此时需要注意的是，利用微扰法反演声速剖面的前提是平均声速剖面下的本征声线和实际声速剖面下的本征声线相近，显然对于海水的声速1 500 m/s，a取100 m/s以内时，平均声速剖面和实际声速剖面的折射率几乎无变化，可以满足本征声线相近的前提。

微扰近似的条件为

$ \begin{array}{*{20}{c}} {t = \int\limits_l {\frac{1}{{\left( {\bar c'\left( z \right) + \Delta c'\left( z \right)} \right)}}{\rm{d}}s} = \int\limits_l {\frac{{\bar c'\left( z \right) + \Delta c'\left( z \right)}}{{{{\bar c}^{'2}}\left( z \right) + \Delta {c^{'2}}\left( z \right)}}{\rm{d}}s} \approx }\\ {\int\limits_l {\frac{{\bar c'\left( z \right) + \Delta c'\left( z \right)}}{{{{\bar c}^{'2}}\left( z \right)}}{\rm{d}}s} = \int\limits_l {\frac{1}{{\bar c'\left( z \right)}}{\rm{d}}s} - \int\limits_l {\frac{{\bar c'\left( z \right) - \Delta c'\left( z \right)}}{{\Delta {c^{'2}}\left( z \right)}}{\rm{d}}s} \approx }\\ {\int\limits_l {\frac{1}{{\bar c'\left( z \right)}}{\rm{d}}s} - \int\limits_l {\frac{{\sum\limits_{k = 1}^{n{\rm{Vectors}}} {{{\alpha '}_k}{{\mathit{\boldsymbol{f'}}}_k}\left( z \right)} }}{{{{\bar c}^{'2}}\left( z \right)}}{\rm{d}}s} = \bar t' - \Delta t'} \end{array} $

(14)

式中：Δc′是一小量。由式 (14) 发现，微扰近似中舍去了传播时间计算式中的Δc′²(z) 项。若Δc′²(z) 比Δc²(z) 大，略去Δc′²(z) 则使得分母变大，因此反演出的结果中分子项c′(z)-Δc′(z) 也将增大，导致反演结果c′(z)+Δc′(z) 减小。

为使反演结果更精确，考虑修正微扰近似引起的误差，可知Δc′(z) 为常数a附近的扰动量。设x(z) 为反演结果的偏移量，令反演得到的Δc′(z) 为在a+x(z) 附近的扰动量，即Δc′(z)≈a+x(z)，代入微扰近似的等式：

$ \frac{{\bar c'\left( z \right) + \Delta c'\left( z \right)}}{{{{\bar c}^{'2}}\left( z \right) + \Delta {c^{'2}}\left( z \right)}} = \frac{{\bar c'\left( z \right) - \Delta c'\left( z \right)}}{{\Delta {c^{'2}}\left( z \right)}} $

(15)

可以得到偏移量满足的等式：

$ \frac{{\bar c'\left( z \right) - a}}{{{{\bar c}^{'2}}\left( z \right) - {a^2}}} = \frac{{\bar c'\left( z \right) - a - x\left( z \right)}}{{\Delta {c^{'2}}\left( z \right)}} $

(16)

解得偏移量：

$ x\left( z \right) = \bar c'\left( z \right) - a - \frac{{{{\bar c}^{'2}}\left( z \right)\left( {\bar c'\left( z \right) - a} \right)}}{{{{\bar c}^{'2}}\left( z \right) - {a^2}}} $

(17)

2 改进微扰法的仿真验证

为了验证本文所提出方法的可行性，做如下的数值计算。数值算例中，假设海底水平，水深50 m，声信号频率10 kHz，声源深度10 m，利用水平阵不同阵元接收信号的传播时间反演海水中的声速剖面。水平接收阵深度40 m，水平阵Ⅰ的阵元距离声源的水平距离分别为13、13.4、13.8、14.2、14.6、15 km；(水平阵Ⅱ的阵元距离声源的水平距离分别为53、53.4、53.8、54.2、54.6、55 km；水平阵Ⅲ的阵元距离声源的水平距离分别为13、14、15、16、17、18 km)，按从近到远的顺序用1~6号对阵元命名，仿真环境的示意图见图 1。

仿真计算中，将真实的平均声速剖面下调20 m/s得到调整后的平均声速剖面，分别使用调整前后的平均声速剖面和对应的经验正交函数进行反演计算，表 1给出了平均声速剖面调整前后的平均声速与实际声速下的声线传播时间差。可以看出，声速剖面调整后的传播时间扰动量的量级为10^-1 s，远大于声速剖面调整前10^-3 s，此时传播时间扰动抗传播时间测量误差干扰的能力得到增强，反演结果的精度得到提高。值得注意的是本文提出的方法相对于传统的微扰法仅仅是在最初的平均声速剖面表示上给定一水平偏移量，以此来改变经验正交函数表示。因此相对于传统的微扰法没有加大计算量，没有牺牲计算时间。

考虑传播时间测量误差对反演算法的影响，假设测时误差服从高斯分布，均值为零，方差为10^-4 s。进行100次声速剖面反演，取反演声速的平均值作为最终反演结果。图 2给出了平均声速剖面调整前后，利用水平阵1各阵元接收信号的传播时间进行反演获得的声速剖面。比较图 2(a)、(b)中5条样本声速的反演结果可以看出，平均声速调整后，反演获得的每条声速剖面均与待反演声速剖面更接近，证明了改进声速剖面反演算法的可行性。然而从图 2(a)、(b)的反演误差曲线可以看出，图 2(b)中的反演结果存在0.3 m/s的水平偏移，这验证了上节中的理论分析。应用式 (17) 对该结果进行补偿，图 2(c)给出了补偿后的声速剖面。可以发现，图 2(c)的反演结果 (实线) 与实际待反演的声速剖面 (虚线) 几乎完全重合。从表 2给出的反演误差可以看出，本文提出的反演算法相比于传统的微扰法有更高的反演精度，验证了此方法的可行性。

为了验证本文提出方法的有效性，下文利用水平阵Ⅱ和水平阵Ⅲ，采用传统微扰法进行仿真研究，反演结果见图 3。且水平阵Ⅱ反演声速剖面的最大误差为1.173 m/s，均方根误差为0.324 m/s；水平阵Ⅲ反演声速剖面的最大误差为0.578 m/s，均方根误差为0.298 m/s。水平阵Ⅱ相对于水平阵Ⅰ而言，阵元间距相同，但距离声源的水平距离更远。与表 2数据比较发现，当阵元间距一定的情况下，增大基阵距离声源的水平距离可以提高反演结果的精度，这是因为声速扰动量敏感性增大的结果。水平阵Ⅲ相对于水平阵Ⅰ有着更大的阵元间距，比较发现，增大阵元间距亦可以提高反演结果的精度，这是因为声波到达各阵元的时间差得到增大，声线轨迹近似导致声线传播时间误差的影响得到降低。

虽然通过增大基阵与源的水平距离或扩大水平阵阵元间距可以降低反演结果的误差，但是传统微扰法的反演精度仍比本文提出的改进算法低。从图 3还可看出，利用水平阵Ⅱ和水平阵Ⅲ反演获得的声速剖面与实际声速剖面相差较大，而本文提出的方法的反演结果如图 2(c)所示，任意一条反演得到的声速剖面都与实际声速剖面极为吻合，验证了本文提出方法的反演性能。

3 结论

本文在传统微扰法反演声速剖面的基础上，通过调整平均声速剖面以增大声速扰动量，增大了声速扰动引起传播时间的扰动量，提出了声速剖面反演改进算法，数值仿真得到以下结论：

1) 声速剖面调整不仅提高了反演算法抗传播时间测量误差的能力，同时也增强了算法对声速扰动的敏感性，有效提高了反演结果的精度。

2) 通过对声速剖面调整引入的声速偏移进行修正，进一步降低了改进算法的反演误差，进一步提高了反演结果的精度。

3) 本文提出的方法相比于传统的微扰法在保证计算时间的情况下提高了反演结果的精度，降低了声波传播时间测量误差对声速剖面反演算法的影响，验证了本文提出方法的可行性。

此外本文提出的改进算法可应用于利用经验正交函数进行反演声速剖面或海洋流速中，降低测时误差的影响。本文仅对此算法进行了仿真分析，海上实验的验证将是我们下一步研究的重点内容。